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ABSTRACT OF THE DISSERTATION 


The Accuracy of Parameter Estimation in 
System Identification of Noisy Aircraft Load 

Measurement 

by 

Jeffrey Kong 

Doctor of Philosophy in Electrical Engineering 
University of California, Los Angeles, 1994 
Professor Kung Yao, Chair 

This thesis focuses on the subject of the accuracy of parameter estimation 
and system identification techniques. Motivated by a complicated load measure- 
ment from NASA Dryden Flight Research Center, advanced system identification 
techniques are needed. The objective of the problem is to accurately predict the 
load experienced by the aircraft wing structure during flight determined from a 
set of calibrated load and gage response relationship. We can then model the 
problem as a black box input-output system identification from which the sys- 
tem parameter has to be estimated. Traditional LS techniques and the issues of 
noisy data and model accuracy are addressed. A statistical bound reflecting the 
change in residual is derived in order to understand the effects of the perturba- 
tions on the data. Due to the intrinsic nature of the LS problem, LS solution 
faces the dilemma of the trade off between model accuracy and noise sensitiv- 
ity. A method of relating the two conflicting performance indices is presented, 



thus allowing us to improve the noise sensitivity while at the same time confin- 
ing the degradation of the model accuracy. SVD techniques for data reduction 
are studied and the equivalence of the Correspondence Analysis (C A) and Total 
Least Squares Criteria are proved. We also looked at nonlinear LS problems with 
NASA F-lll data set as an example. Conventional methods are neither easily 
applicable nor suitable for the specific load problem since the exact model of the 
system is unknown. Neural Network(NN) does not require prior information on 
the model of the system. This robustness motivated us to apply the NN tech- 
niques on our load problem. Simulation results for the NN methods used in both 
the single load and the ’’warning signal” problems are both useful and encourag- 
ing. The performance of the NN(for single load estimate) is better than the LS 
approach, whereas no conventional approach was tried for the ’’warning signals” 
problem. The NN design methodology is also presented. The use of SVD, CA 
and Collinearity Index methods are used to reduce the number of neurons in a 
layer. 


xm 



Chapter 1 
Introduction 


Parameter estimation and system identification are general analytical techniques 
applicable to many engineering problems. The problem can generally be viewed 
as a black box system identification process from which input and output relation- 
ship is identified and determined. Once the system parameters are determined, it 
can be used to estimated the output when a new set of input data are presented 
to the system. As in most cases such as spectral estimation in signal processing, 
statistical analysis and data reduction techniques, our problem reduces to finding 
the solution of the linear system of equations represented by Ax & B, where A is 
the input matrix and B is the output matrix. The parameter x relating A and B 
is then determined using linear algebraic techniques. In the past, normal equa- 
tions approach is mostly used until the development of newer techniques such 
as Singular Value Decomposition, QR factorization, etc. However, there are still 
drawbacks in these methods which we need to understand and analyze. There 
are factors affecting the accuracy of the system identification process which we 
need to address and analyze. Problems that appeared in this process include the 
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redundancies in the data matrices, the need for subset data selection, the noise 
in the data matrices, the selection of system parameters and the accuracy of the 
linear model assumption. Some of these concerns are studied in [30] [9] such as 
using S VD to solve the problem of data redundancies. In [3] and [9] , Correspon- 
dence Analysis is used in data reduction. In [7] [10] [11], the presence of noise in 
both the data matrix A and output matrix B are considered using Total Least 
Squares. Also, we need to address the intrinsic problem of model accuracy and 
its robustness to noise. 

In practical situations the input-output relationship of the black box system 
cannot be accurately modeled a s a linear systems. As a result linear methods 
may not be adequate for these situations. Techniques such as spline approxima- 
tion and nonlinear optimization are widely used in solving nonlinear least square 
problems. Spline approximation is relatively simple in dealing with only one or 
two dimension problem but becomes increasingly complex in multivarite prob- 
lems. On the other hand, the use of nonlinear optimization techniques is not 
always straightforward and simple. Often, the assumption of the type of nonlin- 
ear function is required for the model. Recently, much attention is given to the 
use of Neural Network approach in applications arising from different areas of 
interest. From pattern recognition, function approximation, linear and nonlinear 
programming to singular value decomposition, Neural Network are used instead 
of more traditional methods. In some cases when model assumption is not easily 
obtainable, Neural Network will be most suitable since no model assumption is 
required. This advantage of NN becomes more apparent when nonlinearities is 
present in the data matrix A. Linear methods are no longer adequate in these 
situations, as a results nonlinear least square techniques must be used. Since 
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NN is based on search algorithms such as gradient descent algorithm, the com- 
putation effort required for the training of the NN is generally very intensive. In 
addition there is no definite criteria for selecting which type of NN to be used for 
the problem, the* use of NN sometimes becomes heuristic and ad-hoc in nature. 
However, the ease of use, the lack of necessary model assumptions and the ability 
to handle wide variety of problems makes it a relatively popular method among 
various area of interest. 

This thesis is motivated by a real aircraft system identification problem en- 
countered at NASA Ames-Dryden Flight Research Center. The physical problem 
deals with the wing surface of an aircraft which is constantly experiencing dif- 
ferent loadings during the flight. The ability to estimate these in-flight loadings 
are essential to the understanding and design processes of the wing structure. 
Strain gages are mounted on different parts of the wing which are sensitive to the 
loads. In order to relate the gages’ outputs to the loadings on the wing surface, 
a pre-flight calibration procedure is performed. The calibration stage is simply 
a procedure to obtain the gages’ outputs when a set of known wing loads are 
applied to the wing structure. In Figure 1.1, we apply a known load L on the 
aircraft structure and obtain the strain gage measurements M. From these data, 
we obtain a set of parameters that characterizes the structure. In Figure 1.2, 
during the flight measurement stage, from the in-flight gage measurements M, 
we can then estimate the in-flight loads L. 

Some major issues in the parameter estimation problem are noises present 
in the calibration? and. the in-flight stages; the accuracy of the model (linear or 
non-linear) used to- establish the relationship between gages and loadings; and 
the number of the available gages. We shall consider some of these matter. In 
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Figure 1 .2: Parameter Estimation Process-Calibration Stage 

chapter 2, some known analytical techniques useful for least-squares solution 
are summarized. In chapter 3, the equivalent approaches on the formulation 
of the load measurement problem are considered. In chapter 4, a bound for 
evaluating the effect of noise on the residual of the LS problem is derived. In 
chapter 5, the statistical properties of noise in both the calibration and in-flight 
stage are discussed. In this chapter the relationship between the noise sensitivity 
and the model accuracy is presented, which leads to the development of a new 
methodology for reducing the flight stage noise variance. It is also shown that 
the reduction of noise variance in the flight stage depends on the amount of 
deterioration of the model accuracy. In Chapter 6, the methodology of Total Least 
Square (TLS) [7], technical results and studies are presented with reference to the 
works by Golub in [7] [10] [11]. We will also discuss the method of Correspondence 
Analysis (CA) studied by Greenacre in [9]. In the same chapter, we will show 
the equivalence of TLS and CA as both methods require the minimization of the 
energy of perturbation error imposed on the data matrix. 
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Chapter 7 introduces the nonlinear least square problems together with a 
practical example by the NASA F-lll gage data. The inadequacy of the lin- 
ear methods leads to the use of nonlinear least square or nonlinear regression 
techniques [2]. We also provide a background of the theory and applications of 
Neural Networks (NN) [20] [17] [12]. Different types of Neural Networks axe also 
discussed and analyzed. One of the most useful type of NN is through the use of 
the Back Propagation training approach [27], Since NN is shown to be useful in 
applications such as pattern recognition, function approximation, nonlinear pro- 
gramming, etc., examples and studies relating to these topics are also presented. 
In Chapter 8 we propose the use of NN in nonlinear least squares and pattern 
recognition approaches in solving the system identification problem represented 
by the NASA F-lll gage data. Other than estimating a single load on the wing 
surface, we propose the estimation of a distributed load. A distributed load es- 
timate consists not only of the information of the single load estimate but also 
the individual loads experienced by different part of the wing. Pattern recog- 
nition is used in an effort of separating different loading conditions. Using NN, 
we showed that we could identify different loading conditions during flight stage 
with reasonalbly amount of accuracy. A simulated examples is also presented. 
In the same chapter, we also proposed a new method to estimate the single load 
estimate using NN. Since NN is capable of solving nonlinear least square prob- 
lems, we could apply the techniques to the NASA F-lll data when conventional 
linear methods failed. Results using NN and linear methods are discussed and 
compared. In the last chapter, the performance of the NN is discussed and a 
design methodogly is presented with simulation results. 
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Chapter 2 

Least-squares Estimation 


Least-squares is a general numerical method that can be used for solving a linear 
system of equations such that the sum of squares of the errors is minimized. 
Suppose we have a system of equation 

Axab, ( 2 . 1 ) 

where A is a known mxn matrix, x is an unknown n x 1 solution vector ,and b is 
a known m x 1 vector. For over determined system of equations, we have m> m 
and for underdetermined system, we have m < n. We also assumed all the data 
and solution are real-valued. The least-squares (LS) problem is equivalent to the 
minimization problem of 

mjn ||Ar — 6|| = ||i4x — 6||, (2.2) 

where x is the LS solution for equation (2.1). Three well known methods [7] for 
solving x are summarized below. 
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2.1 Normal Equation Method 

By expanding (2.2), we have 

1 1 Ax — b\\ 2 = ( Ax — b)'{Ax — b) = x'A'Ax — x'A'b — b' Ax + b'b 

= x'A'Ax - 2x>A'b + b'b. (2.3) 

By taking the gradient of (2.3) with respect to x and set it to zero yields 

2 A' Ax - 2 A'b = 0 

and 

x = (A'A)- l A'b. (2.4) 

Thus the normal equation approach to the least-squares solution is given by (2.4). 
This solution exists only if m > n and A is full rank. However this method does 
not yield a numerically stable solution since it requires the inverse of the matrix 
A' A whose condition number is the square of that of A. Thus this method is not 
useable if A is already nearly singular. 

2.2 QR Factorization Method 

Let A be a full rank mxn matrix that can be factored into the form of A = QR, 
where Q is a, mxn matrix of orthonormal column vectors, and R is a n x n upper 
triangular matrix [7}. Then from (2.1), x can be solved by the equations 


QRx = b 

(2.5) 

Rx = Q'b. 

(2.6) 
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Since R is a non-singular matrix, R 1 exists and (2.6) can be written as 


x = R-'Q'b. (2.7) 

While the formal LS solution x based on the QR method is given by (2.7), since R 
is upper triangular, x in pratice can be obtained by using the back substitution 
method once Q'b is determined. In general, the QR method approach for the 
solution of x is more stable than the normal equation approach. 

2.3 Singular Value Decomposition Method 

Let A be a m x n matrix of rank r < min (m,n). Then the SVD of A can be 
expressed as A = U’ZV’, where U is a m x n matrix with orthogonal columns 
(i.e., U'U = /); V is a n x n orthogonal matrix (i.e., VV = V’V = /); and E is 
a n x n diagonal matrix = diag[Ai, ... ,A n ], with Aj > A 2 > • • • > A r > A r+ i = 
■ ■ • = A n = 0 [7]. Then from (2.1), x can be solved by 

x = U'X + Vb, (2.8) 

where E + = diagfAf 1 , ... ,A“’ , 0, . . . , 0 ]. Thus the SVD solution of (2.8) is 
numerically stable even when A is rank deficient. As a result, the SVD approach 
toward the solution of x in (2.8) allows effective reduced rank of the system of 
equations. 
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Chapter 3) 

Two Equivalent Approaches to 
Load Measurements 


There axe two fundamental operations in the load measurement problem on an 
aircraft structure. In the calibration stage as shown in Fig. 1.1, we perform a 
system identification by measuring the responses of the strain gages mounted on 
different locations of the structure from a series of known applied load at various 
specified load points on the structure during the calibration stage on the ground. 
In the flight measurement stage, by using some characterization of the system 
obtained in the first stage, we can predict the actual equivalent load value and 
location from the gage measurements during the flight. Various known successful 
approaches and results have been reported in the past on the load measurement 
problem [14] - [26]. There are two fundamental and intuitively equally justifiable 
linear approaches arbitrarily denoted as Approach 1 and Approach 2, that are 
applicable to the load measurement problem. In Approach 1 , we model the load 
value matrix L as dependent linearly on the influence coefficient value matrix M 
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measured by the gages. In Approach 2, we model M as dependent linearly on L. 
In general these matrices are rectangular, thus it is not immediately clear that 
these two approaches are equivalent. Historically, all the work in [14] [26] were 
based on that of Approach 1 . Now, we shall show that these two approaches are 
indeed equivalent in most cases, and can be proved by the use of the SVD method 
of Section 2.3. On the other hand, if we only use the more conventional and 
previously used normal equation method ([14] [26] also called the linear regression 
technique) of Section 2.1 , then the limitation of this analytical technique can only 
show the validity of Approach 1 when the number of gages n is less or equal to the 
number of loads m(m < n). Through the use of the concept of ’’minimum energy” 
solution, we can show that both approaches are equivalent. There are several 
theoretical, practical, and computational consequences to these observations. 

At the most basic level of understanding, of course, it is theoretically impor- 
tant to know the equivalency of these two seemingly different approaches that 
yield the desired result. At the practical algorithmic operational level, the inad- 
missibility of having the number of gages n greater than the number of applied 
loads m in the calibration stage in Approach 1 is generally not serious. However, 
there are certain conditions in which we want to consider more gages than the 
number of loads in the calibration stage. Conventional normal equation approach 
(i.e., Approach 1) is not possible since A' A needed in the processing is singular. 

When the data from the gages are quite linearly independent, then there is 
no significant numerical difference between the use of the SVD technique or the 
normal equation technique. However, for highly dependent data, there can be 
significant advantages for the SVD technique. Detailed numerical computations 
based on practical observed gage measurements and load values axe necessary to 
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verify their differences. The crucial point is that in all cases, the SVD approach 
is always computationally more costly as well as numerically more stable. For 
typical dimensions encountered in the load measurement problems, the additional 
computational cost of the SVD approach is not of significant concern, when we 
perform only few LS computations. However, when we perform the LS com- 
putations repeatedly (as we shall see in Section 3.3 under the exhaustive search 
method), then the higher SVD computational cost may be objectionable. 


3.1 Approach 1 - Linear Dependency of Load 
Values on Gage Values 

Load Matrix 

Let L G lR mx3 be a load matrix 


L=[L U L 2 ,L 3 ], (3.1) 

where 

Li = L s = [S u S 2 ,...,S m ]', (3.2) 

L 2 = L b = [Bi , B 2 , . . . , £ m ]', and (3.3) 

L 3 = L t = [Tj, . . (3.4) 

Alternatively, the shear, moment and torque of the i-th element can be ex- 
pressed as 



i — 1, .. . , m, and 

(3.5) 

S%Xi , 

i = 1, . . . , 771 , 

(3.6) 
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where ( Xi,yi ) represents the relative position of the load vector. 

Gage Matrix 

Let M € ]R mxn be the gage matrix which is the response of the n gages to 
the m loads in the calibration process, specifically let 


M = [m.i, . . . ,m. n ] = 



(3.7) 


where each m,i, i = 1 , . . . , n, represents the normalized response of the i-th gages 
to the m loads. Let the n x 3 dependency coeifficient matrix b consists of 


b = [ 6 #i, 6 . 2 , 6 * 3 ] j (3-8) 

Then the linear relationship between the load and the gage matrix can be 
expressed as 


Mb 

(3.9) 

M[b,i,b, 2 ,b, 3 ] or 

(3.10) 

Mbti ,i= 1 , • . . , 3. 

(3.11) 


For i = 1, b.i yields the dependency of L\ = Ls, the shear vector to the linear 
combinations of the influence coefficient vectors m, 1 , . . . ,m. n of M in (3.7). Simi- 
larly, 6*2 yields the dependency of L« 2 = Lb and M and b . 3 yields the dependency 
of L, 3 = Lt and M. 

For the cases of m > n, normal equation and Approach 1 

In the calibration process, the matrix M as well as Li, L 2 , and L 3 are available. 
If we define the pseudo-inverse of M € ]R nxm as a matrix from the normal 


12 




equation point of view, we have 


M + = (M'M)~ 1 M I , (3.12) 

where m > n and all columns of M are linearly independent. Then the LS 
solution b of the linear system L ~ Mb using the normal equation method is 
given by 

b.i = M + L u i= 1,2,3, and (3.13) 

b = M + L. (3.14) 

In the flight measurement stage, we observe one 1 x n dimensional gage mea- 
surement vector M (corresponding to the first row vector of M in (3.7)). From 
(3.9), the predicted 1x3 load vector L is given by 

l = [S,B,f] = Mb= MM+L^MiM’My'M'L (3.15) 

= MM + [L. lt L. 2 ,L. 3 ]. (3.16) 

Then the first element of L yields the predicted shear S, the second element 
yields the predicted moment B, and the third element yields the predicted torque 
T as given by 


S = MM + L. U 

(3.17) 

B = Sy = MM + L, 2 , and 

(3.18) 

T = Sx=MM + L. 3 . 

(3.19) 

From (3.19), we can solve for y and x as 


MM+L.a 

y — — and 

MM+L.i 


. MM+L, 3 

X ~ MM+L.x' 

(3.20) 
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Thus, (3. 19) and (3.20) represent the predicted equivalent net shear, bend- 
ing moment location, and torque location of the applied load that yielded the 
measured gage vector M using the normal equation approach. 

For the cases of m < n, normal equation and Approach 1 

For this case, we have an underdetermined linear system. Since m < n, we 
can solve the linear system Mb = L exactly. However, a unique LS solution 
doesn’t exist for the normal equation method. Therefore a ’’minimum energy” 
solution can be chosen for this case. 

Let / ( b ) be an energy function to be minimized with a constraint that Mb = L, 
we can use Lagranian multiplier method such that 

m 

fk(b. k ) = b>. k b. k + i:\ J (m j .b. k -L. k ) (3.21) 

i = i 

= K k b. k 4- \'{Mb. k - L. k ] = b: k b. k + [b'. k M' - L' %k ) A, (3.22) 

where k = 1, 2, 3 ,A = [A* , . . . , A m ]' is the vector of the lagranian multiplier. Now 
setting the gradient to 0 yields 

V/*(6. fc ) = 2 b. k + M'X = 0 (3.23) 

and 

b. k = -jM'A (3.24) 

A = —2(MM')~ l L, k , k= 1,2,3. (3.25) 

As a result, the solution can be expressed in matrix form as 

b = M\MM')~ l L = M + *L. (3.26) 
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SVD Method for Approach 1 

Now consider the use of the SVD technique via Approach 1. Consider an 
alternate form of the SVD of the matrix M with rank p < min(m, n ) as given by 

M = C7„E„ V„, (3.27) 


where 

Um € IR mxp is a matrix with orthogonal columns, 

V M € IR nxp is a matrix with orthogonal columns, and 

'Em € IR pxp is a diagonal matrix with positive singular values (S.V.) denoted 

by <Ti > (72 >...> (T p > 0. 

Let the pseudo-inverse of the matrix M 6 IR nxm be denoted by M ++ , and 
from the SVD point of view we can express 

M ++ = V m E+ U’ m , (3.28) 

where E^ is a n x m diagonal matrix with the inverse positive singular values 
denoted by \/<Jx < l/<x 2 < ... < 1 /<r p > 0. Thus by using (3.28) in (3.9) , the 
LS solution for b satisfying (3.9) via SVD is 

b = M^L. (3.29) 

We note, (3.29) corresponds to (3.14) and (3.26) in the calibration stage of the 
previously considered normal equation technique. 

Then in the flight measurement stage , we have 

L = Mb — MM ++ L. (3.30) 

It is most interesting to note, that the predicted load vector in (3.30) based on 
the SVD technique has the same form as the predicted load vector in (3.16) based 
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on the normal equation technique. Indeed when m>n (i.e., the number of loads 
is greater or equal to the number of gages) and when the gage measurements are 
quite linearly independent, the pseudo-inverse given by M ++ in (3.29) is equal to 
the pseudo-inverse given by M + in (3.9). Thus, in those cases, either the conven- 
tional normal equation or the SVD methods will yield the same predicted load 
values. Of course, when the measurement values are quite linearly dependent, 
then the SVD approach will be better from the numerical stability point of view. 
As mentioned earlier, when n > m, the normal equation method is not applicable 
for Approach 1 since M + in (3.12) is not defined. Through the use of the concept 
of ’’minimum energy” solution, we showed that normal equation approach can 
be used for cases when n > m. In addition, the results of (3.28)-(3.30) under the 
SVD method for Approach 1 are valid in all cases including n> m. 
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3.2 Approach 2 - Linear Dependency of Gage 


Values on Load Values 

Prom a physical cause and effect point of view, it is reasonable to assume that 
the response of the first gage to the m loads is given by 



’ mn " 


' S\ s\y\ S\X\ 


’Cxi ' 

771*1 — 

• 


: : : 


Cl2 


- Tflrnl - 


- S m J/m SmXfn . 


-Ci3. 


= [L.i, L. 2 , L, 3 ]ci = Lc\. (3.31) 

In (3.31), we are describing the gage measurement as a linear combination of the 
form 

mu = SjCn + SiyiCi2 + s»z»ci3, (3.32) 

which depends linearly on the shear, bending moment, and torque. In general, 
for all n gages, we have 

M = [m,i, . . . ,m. n ] ~ L[ci,. . . ,c„] = LC, (3.33) 

where the 3xn dependency matrix C is denoted by 

C — [C, > ■ ■ • j Cn] . 

In the calibration process, M and L are available as before. In the flight measure- 
ment process, we have an observed M given by (3.33) as 

M = LC. (3.34) 
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For the cases of m < n, normal equation and Approach 2 
Let the psudo-inverse of L be denoted by 

f (L'L)~ X U for m > 3 

L + = < v ' “ (3.35) 

[ L(LL')~ X for m < 3. 

By using (3.35) in (3.33), the LS solution satisfying (3.33) is 

C = L + M = (L'L)~ X L'M. (3.36) 

In the flight stage, we have 

M = LC = LL + M. (3.37) 

Since m <n, (MM') -1 exists and we can multiply both sides of equation (3.37) 
by M' and obtain 


MM' = LL + (MM') 

(3.38) 

MM'(MM')- 1 = LL + 

(3.39) 

MM'(MM')- l L = LISL 

(3.40) 

hx3 


L = M(MM')M'L. 

(3.41) 
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For the case m > n, normal equation and Approach 2 

Using the ’’minimum energy” solution of ( LL + )' for the underdetermined 
system we have in the flight stage from equation (3.36) and (3.37) 


m = id 

(3.42) 

M = (LL + )M 

(3.43) 

M' = M'(LL + y 

(3.44) 

( LL + y = M(M , M) _1 M , 

(3.45) 

LL + = M[(M' M' 

(346) 

L = M(M'M) _1 M'L, 

(347) 


where L + L — /3X3. 

SVD method for Approach 2 

Now, consider solving for C in (3.33) by using the pseudo-inverse of L based 
on the SVD representation of L. Specifically, consider the SVD of L as given by 

L = Ul£. l VI (3-48) 

where 

U L e IR mxm is an orthonormal matrix, 

V L € IR 3x3 is an orthonormal matrix, and 

e IR mx3 is a matrix of the form diagian, (Tl 2 , <7l 3 )- 
In particular, we note 

= h. (3-49) 

By using (3.48) in (3.34), we have 

M = LV L EtU' L M. (3.50) 
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Now, M can be expressed in its SVD form of 

M = U m T, m V' m . (3.51) 

Then (3.50) can be written as 

M = L V L Y.tU' L U m Y, m V' m . (3.52) 

' v 

L+ 

For the cases of m < n, SVD and Approach 2 

For m < n, E& G IR mxm , U M € IR mxm , and V M € IR nxm . Then V M V' M = 
/mxm and UmU’ m = U' M U M = / mxm since U M is orthonormal. 

Then equation (3.52) becomes 


M(V u Xt,U' u ) = 

'LL + 

(3.53) 

IfSL, - 


(3.54) 


M++ 


L = 

MM ++ L. 

(3.55) 


Since equation (3.55) in Approach 2 is the same as that of Approach 1, we showed 
that the two approaches is equivalent in this case. 

For the cases of m > n, SVD and Approach 2 
For m > n we look at the expression of 

M « LC. 

The SVD solution C = Uc^cV' of the above LS problem is 


C = L + M. 


Since C € 1R 3xti and CC + = 1^x3, we can write 


CC£ = L+MC+ 

(3.56) 



C + = M ++ L. 

(3.57) 
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For the in-flight stage we have 


M = LC 

(3.58) 

= ic&i 

(3.59) 

M++L h*3 


l = MM ++ L. 

(3.60) 


Since equation (3.60) in Approach 2 is the same as that of Approach 1 , hence we 
showed that the two approaches yield the same solution via SVD. 


3.3 Limits on the Number of Strain Gages 

Due to limitations of in-flight telemetry channels, the number of gages available 
for the measurement of loads are less than that during the Calibration Stage. 
Thus we need to find the best possible combination of smaller number of gages 
for load estimation. Two heuristic approaches have been used to eliminate gages 
in such a way that the resulting increase in estimation error is not very significant 
126]. 

T Value Method 

For each gage, a T value is calculated as 

Ti = &/«, (361) 

where 6, is the i-th element of b and c t is the standard deviation error associated 
with bi [23]. Then the gage with the smallest T value is eliminated. Repeat the 
process until the desired number of gages is left. 
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Figure 3.1: Residuals vs. number of gages used 


MT Value Method 

Again, for each gage, a MT value is computed as 

MTi = (3.62) 

7«^i 

where [m is the mean absolute response of z-th gage, A i is the correlation coefficient 
from z-th gage to all load conditions, and 7 < is the correlation coefficient from z-th 
gage to all other gages. The elimination procedure based on the MT value is the 
same as that based on the T value. 


Exhaustive Search Method 

This method minimizes the residual by exhaustively computing all the combi- 
nations of the desired number of gages. The combination which results in the 
lowest residual will be chosen as the best gage combination for load estimation. 
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A simulation using NASA HWTSS data is carried out to compare the dif- 
ference in performance for different subset selection methods. In fig. 3.1, we 
observed that although MT and T value methods are heuristic in nature, it pro- 
vides a more accurate solution than the solution yielded by the ’’worst gage” 
choice. Indeed, these methods perform quite close to that of the ’’best gage” 
choice method. Since the choice of which gages to be eliminated depends on the 
integrity of the gage data, the reduced gage choice set can be very sensitive to 
noisy data. Thus a more robust procedure is preferred when data matrix are 
corrupted with noise. Other methods that use collinearity index in trying to 
eliminate redundant gages are also studied in works by [30]. 
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Chapter 4 

A Bound for Noise Perturbed 
Systems 


The performance of the estimation using least squares solution depends on the 
accuracy of the underlying model and the integrity of the data. This chapter 
will illustrate how a corrupted data matrix affects the accuracy of the model as 
reflected by the change in the residual of the LS problem. The effect of noise on 
the measurement matrix M will be shown by deriving a bound for the increase 
in residual as a function of the noise. The following preliminary materials are 
necessary for the derivation of the bound [7]. 

4.1 Distance Between Subspaces 

Consider an over determined system where m > n. In specific we let A 6 IR mxp 
and B e IR mxp be two given matrices. Let the domain of A and B be IR P , then 
let T € IR m and Q € IR m denote the ranges of A and B respectively. Then we 
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can define 


cos 9k = max max u'v = u' k Vk, k=\,...,p (4.1) 

uefveG 

IM|2=|M|2=1 

u'lLi — v'vi — 0, i = 1, . . . , k — 1, 

where uq and Do are null vectors (i.e., for k = 1, the constraints u'uq = i/vo = 0 
vanish). Then the set 6 k, k= 1, ..., q satisfies 0 < 9i < 62 < • ■ • < 9 q < 7t/2. 
Furthermore, u,, u, and 9, are called the principal vectors and the principal angles 
of the subspace pair of T , Q [7]. In particular, 9 V is the maximum angle between 
two basis vectors in subspace T and G • Furthermore, the distance between T 
and G reduces to 

WQjrQ'yr - Q g Q'g || = ~ cos 2 9p = sin 9 p , (4.2) 

where the columns of Qjr and Qg are the orthonormal bases for T and G respec- 
tively. It can be shown [7] that 

COS0* = <Tk(Q'jrQg),k = l,...,p. (4.3) 

where (JkiQ'jrQg ) is the kth largest singular value of the matrix Q'jrQg. 

4.2 A Bound for the Noise Perturbed Resid- 
ual 

Let [ui , ..., Up) be a set of orthonormal basis vectors in IR m that spans the subspace 
M of the strain gage measurements under the Calibration Stage and ,Up] 
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be a set of orthonormal basis vectors in IR m that spans the subspace M of the 
gage measurement with noise. Specifically, the noise E is defined by 

M = M + E. (4.4) 

Furthermore, we denote Qm = [tii, . . . ,Up] and Qm = [vi, . . . ,v p \. Then QmQ'm 
is a projection operator onto M and QmQ'm a projection operator onto M. 

Under ideal (noise-free) condition, the ideal LS solution b associated with the 
measurement matrix M and the calibration load Lcj, is obtained from 

min 1 1 Mb - Lcab \ | 2 = 1 1 Mb - Lcab \ | 2 - (4-5) 

6 

Then L can be defined by 

L = Mb = Q M Q' M L^b- (4.6) 

Similarly, under noisy condition, the LS solution b associated with the noisy 
measurement matrix M and the calibration load L^b is obtained from 

min || Mb - LcobH 2 = || Mb — Lcab\\ 2 - (4-7) 

6 

Similarly, L can be defined by 

L = MQMQ'MLcab • (4-8) 

Thus the difference between L and L becomes 

L — L — ( QmQ'm ~ QmQ'm )^ cab (4-9) 

and its norm squared is bounded by 

|| L — L|| 2 = \\(QmQ'm ~ QmQm)^^ II — ll^'cabll sin 0 P , (4-10) 
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where || (<?„%, - <?*<%) || 2 = sin 2 *, = 1 - a^(Q' u Q a ) from (4.2). 
For simplicity, denote 

a p( Q' mQm ) = Op* . 

Also denote the noise-free and noisy residuals by 


(4.11) 


c Lcab I4 and 

e = L^-L. (4.12) 


Then (4.9) becomes 

L — L = (Lcab — L) — ( — L) = e — e = (QmQ'm ~ (4-13) 

By applying Triangle Inequality, we obtain 

l|f|| 2 <l|f || 2 + l|i«tl| 2 sin 2 « p (4.14) 


or 

||e|| 2 - ||f|| 2 < ||L„i|| 2 sin 2 0 p = ||Z,^|| 2 (1 - <r p * 2 ). (4.15) 

Thus (4.15) yields a bound on the noise perturbed residual. 

The bound in (4.15) can further be expressed in a more useful form if a v * can 
be expressed in terms of the perturbed matrix M and the noise matrix E. As a 
result, the remaining section will be dedicated to this objective. 

Suppose the matrix M and M have QR factorizations as 

M = QmRm (4-16) 

and 

M = QmRm- (4-17) 
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From (4.7) we have M = M + E, after applying QR factorization to M and 
M, we obtain 

QmRm = QmRm + E (4.18) 

Qm = QmRmR m + ERfii 

= Qm + QmRmR'm\I - RmRm) + ER&- 

If we denote 

Ei = QmRmR~m(I ~~ RmRm) 4 - ERjJ, ( 4 - 19 ) 

then (4.19) becomes 

Qm = Qm + E, (4.20) 

Q'uQa = Q' u (Qm + E 1 ) = Q' u Q u + Q' m E 1 = I + Q' u E i . 

Furthermore if we denote 

E 2 =Q' M E 1 , (4.21) 

then (4.21) can be expressed as 

Q'mQm = I + E 2 . (4.22) 

From properties of the theory of perturbed singular value(p.284 of [7]), we have 
\<r k (D + F) - a k {D)\ < a x {F) = \\F\\ 2 , for k = 1, . . . ,p. (4.23) 

Combining (4.21) and (4.23), we have 

Wp(Q'mQm) ~ a p(I)\ = M 1 + E *) - <7 P (J)I < ^ 1 (^ 2 ) = ll^ll- (4.24) 
MQ'mQm)- 1| < M 
Wp*-1\ < mw 

(l-<r P *) < 1 1 E 2 1 1 , for <7 P * < 1 

< 7 p* > 1 - ||£s||, for a p * < 1. 
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From (4.21), we have 


Ek = QhQmRuRtJV - «M«») + Q«ER- a ' (4.25) 

= R u R~ a ' (I - R a R- M ') + Q'mER-J 
= RmR& - I + Q’mER'm. 

Taking the norm squared of and invoking THangular Inequality, we obtain 
|[fy| ! = \\R M R- a '-I + Q' M ER- a '\f (4.26) 

= IlflM^'ll 2 + 1 + IIEfi^’H 2 . 

If we apply QR factorization to M and M in (4.19), we obtain 

M = M-E (4.27) 

QuRm = QmRm ~ R 

QmRmRm = QmRmRm ~ RRm 

\\R M R^\? = II 2 

< ll««<?flll 2 + \\ER- a ' II 2 < 1 + ll^«'ll 2 

since ||«'„ll = 110*11 = 1 and HOmOsII 2 < I IQmII*II0/bI I 2 = 1- 

Now Rm is an upper triangular matrix, by defining /i we have 

M = ll*«'l| 2 (4-28) 

= \\diag(\ln l/r,)l| 2 (4.29) 

= max 1 /rf, (4.30) 

where r, for i = 1, ...,p are the diagonal element of the matrix R&- 
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Then combining (4.27) and (4.28), we have 


ll^ll 2 < \ + \\ER- a '\\*+l + \\ER- a '\\ 2 (4.31) 

= 2(1 + HE/^'ll 2 ) 

< 2(1 +^||E|| 2 ). (4.32) 

(4.33) 

In summary, we have the following results 

|| e -|| 2 -||e-|| 2 < ||i^l| 2 (l - <V» 2 ) (4-34) 

op* > 1 - ||£ 2 ||, for (T p * < 1 (4.35) 

ll&ll < ^(TTiW). (4.36) 

(4.37) 


As a result, the bound can be expressed as 

l|€-|| 2 - ||e|| 2 < ||£~*II 2 {1 - [1 - v^O + Hiil)) 2 }. (4.38) 

Thus, we have derived a bound for the change in residual due to the error matrix 
E. It is interesting to see that the bound depends on the value of fi (which 
is the largest value of the diagonal elements of matrix R-m found from the QR 
factorization of the noisy data matrix M) and the largest singular value of the 
noise matrix E. 

Now since || • H2 || • ||f> equation (4.38) can be expressed as 

||e|| 2 - ll^ll 2 < IM 2 {1 - [1 - v/ 2 (l + l*||B|| 2 p)| 2 }. (4.39) 

Thus we have related the norm of the noisy M residual c to the norm of the 
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noiseless M residual i and the norm of the matrix error E. This bound pro- 
vides an analytical relation between the norm of the noiseless and noisy residu- 
als (corresponding to M) thereby allowing us to estimate the norm of the noise in 
the data matrix when the other two terms in (4.39) are determined or estimated. 
Although it is only an analytical bound whereas the actual behavior of the system 
is not accurately determined, it allows us to have an understanding of how these 
quantities are related. 
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Chapter 5 

Noise Sensitivity Analysis for 
the Flight Stage Load 
Estimation 


5.1 Effects of Errors on Load Prediction in Cal- 
ibration and Flight Stage 

During the calibration process, the gage measurement observed can be repre- 
sented as a matrix of 

M C = M C + E C , (5.1) 

where M c £ IR mxn is the true gage measurements and E c € lR mxn is a matrix 
with its elements as the error associated with the gage measurement. In order 
to estimate the flight load, the estimation parameter derived from the noisy 
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calibration matrix M c can be written as 


b = b + b E = (M c + E c )+L c , (5.2) 

where b E IR n is the noiseless estimation parameter, bs is the estimation param- 
eter error and Lc € IR m is the calibration load vector. During the flight stage, 
the observed noisy output gage measurement can be expressed as a matrix of 

M 0 = M 0 + E 0 , (5.3) 

where E 0 e IR" is a matrix with its elements as the error associated with the 
output gage measurement defined as M 0 e IR". 

Thus the final noisy estimated load taking into consideration of the calibration 
and flight stage can be expressed as 

Mob = (M„ + E 0 )'(b + b E ). (5.4) 

The resulting estimated load can also be expressed as 

L 0 = {M 0 + E 0 )'{b + b E ) = M'jb + M' a b E + E' 0 b + E'J)e 
= L 0 + M'J)e + E'J) + E' 0 bE 

= L a 4- toci (5*5) 

where L 0 is the noiseless load estimate and is the estimate that we wish to obtain. 
Similarly toe is the corrupting error term which reduces the accuracy of the load 
estimate. 

If we examine the variance of toe, we could gain some insight orr how to 
reduce the effect of this noise on the L. The variance of toe ca* 1 be expressed as 

var{toc } = var{MobE + E'J) + EJe}- (5-6) 


33 






Calibration Stage 
(LS methods ) I Calibration 

load 


Noisy prediction 




If we assume that the output error E a and b E are uncorrelated, which is reasonable 
since input and output noise are usually uncorrelated, then 

var{E' 0 b E } = 0. (5.7) 

As a result, (5.6) can be written as 

var{coc} = var{M' Q bE } + var{E' 0 b}. (5.8) 

(5.6) indicated that if we want to minimize the var{eoc}, we can minimize 
the first and second term in (5.8) independently. In other words, we could reduce 
the combined effects by independently reducing the noise variance during the 
calibration and flight stage. 

5.2 Properties of Flight Stage Noise and Load 
Estimate Error 

In the last section, we have shown that the overall noise variance can be reduced 
by independently reducing the noise in the calibration and flight stage. For 
simplicity, we now assume the noise in the calibration is zero( i.e. E c = 0). Then 
we will study the effects of noise on the load estimation accuracy during the flight 
stage. As mentioned before, the calibration load L c and the gage measurement 
M are related by the system of linear equations, 

L c a M c b. (5.9) 

The estimating parameter b — M+ L c will then be used to predict the amount 
of load during the flight stage. The flight stage noise, which can be modeled 
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as an additive noise E 0 = [ej, . . . , e n ]' as in (5.3) on the gage readings M a = 
[mi, . . . ,m n ], is generally present in the measurement system . As a result, the 
noisy flight stage load estimate L a can be expressed as 

L 0 = M 0 b = (M„ + E 0 )b = L a + t D , (5.10) 

A A 

where L 0 is the noiseless load estimate and t a = E a b is the excessive load estimate 
error due to the noise E 0 . 

In order to study the corrupting effects of the load estimate error c 0 , we need 
to examine some of its statistical properties. 

If we assume that the mean of the flight stage noise E a is zero, the expected 
value of t 0 can be expressed as 

E{t a } = E{E' 0 b} = 0. (5.11) 

Then the variance of e 0 is 

Var{t 0 } = E{(t 0 -E{t 0 m = E{tl} = E{(E' 0 b) 2 } 

= E{b'E 0 E' 0 b} = b'E{E 0 E' 0 }b = b'X EE b, (5.12) 



/ 

"ejei ... eien 


where H EE = E{E 0 E' 0 } = E < 


; • . ; 



k 

- 

J 


It is often desirable to reduce the variance of t c since a small variance suggests 
that the corrupting effect of noise on the load estimate will be less severe in the 
flight stage. However the variance of t 0 depends only on the amount of noise in 
the flight stage and the estimator b obtained from the calibration process, it is 
impossible to reduce the variance of e 0 without making changes on the vector 


36 



b. Nevertheless, the accuracy of the predicted load will also be affected when 
changes are made on b. In other words, the reduction of the noise sensitivity of 
the estimate is at the expense of the model accuracy. Therefore, it is important 
to optimize the amount of tradeoff between noise sensitivity and model accuracy 
when noise sensitivity reduction is absolutely necessary. The following sections 
will examine two possible solutions to this problem. 


5.3 The Sensitivity Measure for Noise Perturbed 
Systems 


Since the variance of e 0 is an indicator for the intensity of the noise, it can be 
interpreted as an intensity index of noise sensitivity. If white gaussian noise is a 
reasonable assumption(i.e var{e Q } = b'b<r* a = \\b\\ 2 ^t), the sensitivity index can 
be defined as 


'to 



(5.13) 


where o\ is the variance of the flight stage noise E 0 and a\ o is the variance of to- 


5.4 The Trade-off between Model Accuracy and 
Noise Sensitivity 

In order to examine the effects of the reduction in noise sensitivity on the model 
accuracy, we need to look at the noisy load estimate from (5.10) such as 

L 0 = Mob = {M 0 + E 0 )'b = L 0 + e a . (5.14) 
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5.5 Single Coeifficient Noise Sensitivity Reduc- 


tion(SCNR) 

In theory, we would like to have a small valued noise sensitivity index. As a result, 
the reduction of S to = ||6|| 2 in (5.13) requires that the value of the norm of the 
estimator vector b = [Si, . . . , 6„]' to be reduced. It can be achieved by changing 
the values of one or more elements of the vector b. The SCNR method discussed 
in this section only changes one element of bi such that the new estimator vector 
becomes 

A A A A 

bscNR = [feu • ■ • i (kbi , . . . ,6„] , (5.15) 

where c* < 1. Then the decrease of the noise sensitivity AS can be defined as 

AS = ||S || 2 — HSscnrII 2 = (1 — c 2 )S 2 > 0. (5.16) 

Now, we define the noiseless new load estimate(when the output stage noise 
E a = 0) as 

L 0 'SCnr — M'J)scnr = TUibi + . . . + rriiCibi + . . . + Tnnb n (5-17) 

and the new load estimate obtained from the noisy measurement rh 0 as 

L 0 ,scnr = M Q bscNR — + . . . + fhiCibi + . . . + fhnb n . (5.18) 

The change in load prediction accuracy A A with respect to the noiseless load 
estimate L a can then be expressed as 

A A = L 0 — L 0i scnr = (1 — Ci)biTrii. (5.19) 

Combining (5.16) and (5.19), we have 

A A = (bi - yjPi - AS)m, 
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and 


AS = Sf - - bi) 2 = AA(2± - ^). 


771 , 


77li 


(5.20) 


From (5.20), if we fix the change in model accuracy A /l, the decrease in the noises 


sensitivity AS can be maximized by scaling the i-th gage prediction vector b to- 
dbi such that % - f- bi) 2 in (5.20) is the largest among gages i = 1, . . . ,n- 
The maximum AS for a fixed AA = 6 A attained can then be expressed as 


ASscnr — max ~ M 2 }- (5-21) 

o<i<n l. rrii J 

AA = 6A 


Figure 5.2 below shows a plot of the normalized ASscnr versus AA for thff 

NASA Wing data, where the definitions of 

ASscnr 


AS 


‘Ascnr — 


A A = 


lift I 


and 


AA 


(5.22) 

(5.23) 


It indicates the amount of improvement in noise sensitivity ASscnr at a given 
level of accuracy deteriorating A A. 


5.5.1 Sub-optimal Single Coefficient Noise Reduction for 
practical applications 

The SCNR described in the previous section requires the knowledge of rru in 
order to maximize AS with respect to all gages. However, only tti* = tth 4* e, ia 
observable at the flight stage, the exact value of m, is not known. As a result, 
we can only use rhi in maximizing AS , such that (5.22) becomes 

ASscnr = max AS = max b 2 — ( bi) 2 . (524) 

&cyv " 0<i<n 0 <i<n * TUi 

AA = 6 A 


39 





Since to* = rrii + => rrii = fhi — e<* , (5.25) becomes 


(5.25) 

A possible sub-optimal method can be used by restricting the amount of noise Ci 
to the r<7 eo criteria. The max-min criteria can then be expressed as 

maxf min {AS\ra eo criteria}]. (5.26) 

0<*<n l -r<7 eo <e oj <ra eo l ' e ° n 

The above equation (5.26) can be interpreted as the maximum of the minimum 
of the A S' under the rcr ee criteria. By setting = ra ee , the max-min condition 
in (5.26) becomes 

A. Sscnr-so — max AS = 

0<i<n 


Thus we have proposed a sub-optimal method for the SCNR procedure(SCNR- 
SO), which will not perform as good as the theoretical SCNR method. However, 
the analysis of SCNR and the subsequent SCNR-SO methods gives us an under- 
standing of how the model accuracy and the noise sensitivity are related. It also 
provides an analytical bound for the region of operation in which SCNR-SO will 
perform better the original method. 

5.5.2 Performance Analysis of the SCNR procedure 

In order to have a better understanding of the characteristic of the new proce- 
dure, it is necessary to compare the performance of the proposed procedure and 
the original method of load prediction. As stated before, the presence of normal 


£2 / £ ^2 

max 6, — (— bi) 

0<t<n ra i — ra eo 




AA 


mi * -ra eo 


- bi) 2 . 


(5.27) 


. _ C'n , LA/i £ X o 

AS scnr = max AS = max o, — (— bi) . 

0<»<n 0<»<n * 'to* - eoi 

AA = 8 A 
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Figure 5.3: Probabilty Distribution for load estimate L a 

gaussian noise at the flight stage is assumed. One way of comparing their per- 
formance is by evaluating the probability that the noisy estimate L 0 lies within 
a specific range for the original and the new SCNR method. Using general sta- 
tistical methods for establishing confidence interval in normally distributed noise 
model, we can use the interval of ±ra 0 . (For example, a 99.73 % confidence inter- 
val corresponds to r=3) Figure 5.3 shows the areas for evaluating the probabilities 
for the two methods. 
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Now, the probability that the noisy estimate L a lies within the range using 


the original method is 


» - 1 ft'o+r<T t0 X* 

P{original} = P{L 0 - ra to < L 0 < L 0 + ra to } = — / exp(-— )dx, 

Z7T JL 0 -r<r t0 Z 

(5.28) 

where <r to is the standard deviation of the flight stage load error e c in (5.10) and r 
is an integer depending on the level of confidence. Similarly, the probability that 
the new SCNR method noisy estimate L 0> scnr lies within the same range is 

_ ^ 1 f r ** 

P{SCNR} = P{L 0 - ra to < L 0<S cnr < L c + ra to } = — exp(-y)dx, 

(5.29) 


where 


-7-(7 t< , - A A 

<Tio,SCNR 


A A 

ll^scwnll IlfcsCNHlke,, 


(5.30) 


r 2 * = 


r<7 fo - A A 

<7t ol SCNR 


= -jm (5.3i) 

||^SC)V/i|| II^SCN/ilkeo 

Then if the condition P{SCNR} > P{original } in (5.28) and (5.29) is satisfied, 
the probability that the SCNR noisy estimate lies within the specified range is 
higher than that of the original method. In these cases, the new method is much 
superior than the old method. In other words, the noisy estimate is more likely to 
be confined in the specified range when using the new method. Figure 5.4 below, 
which is generated from NASA data, indicates that the SCNR method performs 
better in a significantly large region of operation and especially when large flight 
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stage noise variance cr £o is present. From (5.27) and (5.28), the probability for 
the SCNR-sub-optimal method can be written as 


where 


and 


(5.32) 

(5.33) 

(5.34) 
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Figure 5.4: Probabilty vs AA for NASA data 

In summary, the new procedure essentially generates a biased load estimate 
with a reduced noise variance whereas the original method gives a non-biased 
estimate and a larger variance. Since both of the two methods have desirable 
characteristic in the load estimation, one has to decide how to make trade-off 
between the noise sensitivity (variance) and the model accuracy (bias). Generally 
if a large amount of noise is expected, the new SCNR could be used to reduce 
the corrupting effects of the flight stage noise by confining the deviation of the 
load estimate. 
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5.6 Generalized Method for Reduction of Noise 
Sensitivity (GNR) 

In (5.15), the variance of the error e Q is reduced by changing 6, to Cib t . Now we 
will generalize the method by using the scaled vector be nr such that 

be nr = Be, (5.35) 

where B = diag(b u ...,b n ) and c = [ci , . . . , c*]', c* < 1. Then the new load 
estimate using GNR is 

L 0 ,gnr = M'J)gnr = M' 0 Bc. (5.36) 

From (5.16) and using (5.35), the decrease in noise sensitivity is 

AS = ||6|| 2 — !|£>gw,r| | 2 = (1 — c)'BB( 1 — c). (5.37) 

Also from (5.19) and (5.36), the change in the noiseless load prediction accuracy 
can be written as 

A A = L 0 - L 0i gnr = M’ 0 B- 1 - KBc = Mo'B( 1 - c). (5.38) 

In order to reduce the noise sensitivity by using the GNR method, we need 
to maximize the change in noise sensitivity AS with respect to a fixed A A such 
that 

max AS = min d BBc. 

C C 

M' 0 B{\ - c) = AA (5.39) 

(5.39) can be solved by the Lagrangian method with 

G(c) = c'BBc + \[M' 0 B{\ - c) - A A]. (5.40) 
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Taking the gradient G(c) with respect to c and set it to zero gives 


VG(c*) = 2 BBc * +\'BM 0 (-1) = 0 


(5.41) 


and 

c* = -\b- 2 BM q = - \B~ l M 0 . (5.42) 

/ Zi 

Putting c* into the constraint equation gives 


M' 0 B{ 1 - c) = Ad, 
M’ 0 B{\ + ^B~'M 0 )=&A, 
■ .(AA-MjBt) n (AA-L) 


(5.43) 


Then c* can be expressed as 


X° Ad. 

c * =( mr )B M ° ■ 

Therefore the change in the noise sensitivity A Sgnr can be written as 


(5.44) 


1 


A ScNR — ~ 2 ~( a to - a t 0 ,GNfd — 


e 0 


(L 0 - Ad) 2 

\m \ 2 


(5.45) 


5.6.1 Sub-optimal Generalized method for Noise Reduc- 
tion 

As stated in the previous section, the value of M a in (5.45) is not known. There- 
fore, we need to tolerate a sub-optimal procedure in practical applications. The 
change in noise sensitivity AS using GNR is 

A Sgnr = ll&ll 2 - (L °~ ^ ^ ■ (5.46) 
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Similarly from (5.45), the sub-optimal GNR using the ’’Maximum-minimum Cri- 
teria” gives 

AS G nr-so — max min AS. (5-47) 

C — rjr eo <£o<r<r eo 

Applying the rcr eo Criteria, (5.47) becomes 

(L 0 + ra eo \\b\\ - AA) 2 


ASgnr-so — 


\\M 0 \\ 2 - nr 2 al 


(5.48) 


Co 


5.6.2 Performance Analysis of the GNR procedure 

The probability that the GNR method noisy estimate L 0> gnr h es within the ra eo 
interval can be expressed as 

A J rr2* 

P{GNR-SO } = P{L 0 -n r to < L 0 < L 0 + n r t J = — j exp-— dx, (5.49) 


where from (5.30) 


and 


r i,2* — T 


m 


A A 


t>GNR ~ 


II^GNilll ||&GWrt|ke 
(. L 0 - A A) 


\m 


(5.50) 


(5.51) 


Similarily for suboptimal GNR method, the interval for the evaluation of the 
probability in (5.49) can be expressed as 

||6|| AA 


r 1,2 * * = T 


||5gNR s oII II^GNfi-Solk, 


e 0 


where 




{Lo + raeMzMl 

\\M 0 \\ 2 ~ nr 2 a 2 o 


(5.52) 


(5.53) 
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5.7 Conclusion and Remarks 


This section provided a study of how the trade-off between variance(or noise 
sensitivity A S) and bias(or accuracy A4). We also presented the new methods 
for an effective reduction of noise variance when a given amount of accuracy 
deterioration A .4 is allowed. This problem is also studied in a different setting in 
the field of statistical analysis. In Ridge Regression [7], the problem of constrained 
LS is solved. It solves the LS problem of the min x || Ax — b \ | with the constraint 
that | |x| | < a. However, in this problem the idea relating the prediction accuracy 
and sensitivity is not explicitly used. Thus our study provided an alternate 
approach for analyzing how variance and bias are related to each other. 
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Chapter 6 

Total Least Squares and 
Correspondence Analysis 


The Total Least Squares(TLS) method was popular in signal processing appli- 
cations recently, while the correspondence analysis technique was created over 
twenty years ago and used in applied statistical and data analysis. These two 
seemingly distinct topics developed independently among their own researchers. 
In this section, we provided the basic reformulation and analytical and geomet- 
rical tools to prove the equivalency of these two basic and useful methods. 

In this section, we consider the total least squares estimation method well 
known in numerical analysis and modem signal processing as well as the corre- 
spondence analysis technique encountered in applied data and clustering analysis. 
Due to historical reasons of development, each of these two subjects have gen- 
erally been formulated with its own notations and solutions. However, upon 
more detailed consideration, both optimization problems reduce to the appli- 
cation of SVD technique for the respective solution. Indeed, upon appropriate 
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pre-processing operation of translating the centroids of the data matrices to the 
origin, these two methods are equivalent. 


6.0. 1 Total Least Squares(TLS) 

Least squares techniques based on the statistical assumption that only error is 
present in the vector to be estimated. Total least squares, which is introduced by 
Golub and Van Loan, attempts to generalize the solutions by allowing both ob- 
servation matrix and estimator vector to be corrupted by errors. This assumption 
allows us to have a more consistent solution since errors are usually presented in 
all data. The problem of TLS can be stated as follows. For an overdetermined 
system of equations of the form Ax B, errors occured in both A and B. The 
TLS solutions for x solves the perturbed problem of 

(A + AA) = B + A B, 

with the constraint that the F norm ||AA, Af?]||f is minimized. It can be seen 
that the classical LS problem solves the same problem above except that no error 
is assumed in the matrix A(or AA = 0). 

6.0. 2 Correspondence Analysis(CA) 

In statistical and data analysis, there are myraids of analytical, graphical, and 
intuitive methods for performing data reduction, clustering and display of statis- 
tical properties. The CA techniques was proposed by Benzecre [3] and studied by 
Lebart [19], Greenacre [9] and others. This technique was clearly motivated by 
several standard statistical multivarite techniques of principal component analy- 
sis and discriminant analysis. It provides an analytical method of displaying of 
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a large centered matrix in low dimensional spaces such as the two dimensional 
space. 

Graphical and Numerical Examples 

For a given data matrix X E IR mxn , a basic problem of data reduction is to 
determine and eliminate the set of redundant columns(or rows). If two vectors 
are strictly independent, one of the vectors can be eliminated. However, the 
presence of noise will make the redundant vector strictly independent. The con- 
cept of collinearity in [25] quantify the amount of linear dependencies for a given 
column (or row) vector relative to all other columns (or row) vectors of X. The 
collinearity index of matrix X can be defined as 

Kj = IIZ^IIM, tt roj, 3 = 1, • • • (6' 1 ) 

where X + is the pseudo-inverse of X. 
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We now look at a physical load problem when the 18 x 12 strain gage mea- 
surement matrix X is given by X = 
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These data are obtained from a NASA/Ames hypersonic wing test structure 
(HWTSS) load measurement experiment. The columns represent the output 
responses of 12 strain gages located at the wing and fuselage region of the aircraft 
when 18 known input load conditions are applied on different locations of the wing 
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Figure 6.1: Correspondence Analysis of NASA HWTSS load data 


structure in the calibration stage. In practice, the number of gages available 
during a flight is usually much less than that available during the calibration. 
Thus, from the data in matrix X , we are motivated to determine the redundant 
gages that can be eliminated. Evaluations Collinearity indices are computed as 
k = (24, 4.3, 52, 7.8, 34, 6.2, 29, 7.5, 32, 8.3, 32, 9.5). Applying the CA 
technique gives S = diag(.61, .31, .18, .18, .083, .067, .051, .031, .011, 6.6- 
10 -3 , 5.6 • 10 -3 , 7.2 • 10 -17 ). The coefficients of expansion of the 12 column 
vectors for the two dominant singular vectors are shown in Fig. 6.1. 

From the collinearity index k, we observed that all the odd numbered column 
data vectors are fairly collinear. From Fig. 6.1, we see that these vectors are 
grouped closely together. However, the even-numbered column data vectors are 
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Figure 6.4: Mean translated data vectors 


is zero. In general, the above condition need not be satisfied for arbitrarily P\ 
and P 2 . Now, consider the column zero mean translated version of A and 6. 


>1 = 


(ai + q 2 )/2 

h — 

(6i + 62)/2 

. (ai + a 2 )/2 

i 0 — 

.(61+62)72. 


(6.5) 



'(ai - a 2 )/ 2' 


and A = A — A = 

(ai - a 2 )/ 2 

,6 = 6-6 = 



'(ai - a 2 )/2" 


Pi = 

. (6i - 6 2 )/2 . 

,A = 


(by - fe)/2 

(a 2 - ai)/2' 

( 62 - 60 / 2 . 


( 6 . 6 ) 


After mean translation, P\ and P 2 are anti-symmetric. That is, Pi(l) = 
— P 2 (l) and Pi (2) = P 2 ( 2). Clearly, the optimum straight line runs from P\ to 
P 2 , with Z\ + Z\ = 0 and d, = 2||Pi|| 2 = 2P 2 2 - For any other alternate straight 
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line, Zf 2 + Z '2 > 0, and 


d! 2 = (d'/2 + <f/2) = sj{d/2f - Z[ 2 + yj (d/2) 2 - Z ; f < d/2 + d/2 = d. (6.7) 
Thus, any alternate straight line yields inferior TLS and CA results. 

6.1.2 TLS and CA Criterion for general dimensions of 
data matrix 

Correspondence analysis can be used to represent data collected in matrix form 
in a more compact format by projecting the raw data onto a lower dimensional 
subspace. 

Let A = (dij) e IR mxn with m > n the original data matrix. For simplifi- 
cation purposes, which wil become clear later, let us assume that ]C£li = 0 
(i.e. A is ’’centered”) and let the expression of A = UTV' be its singular value 
decomposition. 

Now, let Q = G IR”* -7 be a matrix with J < n orthonormal 

columns and H = AQ = (hi,), with hij = Oi.q^\ the matrix with its elements as 
the the projections of A onto Q. 

In this context, Q represents the basis-matrix of a J-dimensional subspace. 
For the purposes of correspondence analysis, it is reasonable to vary Q so that 
the projection coefficients are as separated as possible, with respect to each 
basis vector q^\ 

In order to satisfy this requirement, one can perform the maximization of the 
following quantity over the space spanned by the columns of Q such that 

J m m 

d = 

j = l t=i k=\ 
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J m m 

= EEM-2 hiihi + hl,) 

j=l i= 1 k—1 

J m J m m m m 

j = \ t=l j=li=lfc=l j=l fc=l 

J TTI J TTI TTI 

= 2rn^^h* j -2^Y^h ij Y,h kj 

j=l i=l j=l t=l i= 1 

mm m Jj2 

= 2">EE«rEi) 

j=l i=l fc=l m 

J m 

= 2m ^ ^2(hij — Hj) 2 , 
j= i *=i 

where //j = 52™! hij/m. Since A is ’’centered”, it can be shown that Hj = 0, V). 
Therefore 

J m 

max D = 2m max ^ ^ A?i (68) 

j=ii=i 

= 2mmax ||/f|||., (6.9) 

where || • ||f is the Frobenius norm of a matrix. 

The above maximization can be solved by the use of the singular value de- 
composition of A. In fact, 

mil = \\aq\\1 = \\uy,v'q\\ 2 f 
= \mv% 

= < Y, a i- 

i=i 

The equality is obtained when Q = Vj = [u. i, . . . , v.j ] , hence 

V = Q'V = [Ij,0j x{n .j ) ] (6.10) 

and H = AQ = AVj = UUV'Vj = [<7iu.i, . . . ,<jju.j\. This results are known in 
the area of correspondence analysis. 
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Now, let [Q, Q L ] G IR nxn be an orthonormal expansion of Q. Then we have 
11*1! = 

j = 1 

= j^q^A'AqW- £ 

i=l j=J+l 

= £« w v'SV«°>-i|.4q- l ii! 

j=l 

= £X>to 0, V*) J - 11^11! 

j=l fc=l 

= i>mwx. 

k=l 

where E = diag(a x <r n ) G IR nxn . As a result, the maximization of D is 

equivalent to the minimization of 

min||AQ J -||^ = min \\A(I n - QQ')\\ 2 F = min || A - A\\ F = min ||AA||^, (6.11) 
where A - HQ 1 and A A = A — A. 

The maximization problem stated above(or minimizing D) is therefore min- 
imizing the energy of the perturbation imposed on the original data matrix A. 
This perturbation reduces the rank of A to rank(A) < J when the vector X in 

a 

the orthogonal subspace spanned by Q ± solves the equation AX — 0. 

Case for Hj ^ 0 

Now consider the case when A is not ’’centered”, Hj ^ 0 or 7^ 0 

(which is a more general case), then we can write 

max D = 2m max \\H — —l m H\\ 2 F 
q m 

= 2m max ||(/ m - —l m )H\\ 2 F 

TYl 

= 2mmax\\H 0 \\ 2 F 
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= 2m max \\PH\\ 2 F 
= 2m max \\PAQ\\ 2 f , 

where l m is a square matrix of all ones, P — I m — ^l m and H 0 = PH. Thus we 
can still carry out the same data reduction by using A 0 = PA, a ’’centered” form 
of A. 

Case for shifted H criteria 

In order to justify the equivalence of TLS and CA, we need to look at a 
shifted H criteria. Now we consider the shifted projection of A onto the subspace 
spanned by columns of Q such that 

H = AQ + B. (6.12) 


Instead of maximizing D with H — AQ as before, we maximize the shifted 
H = AQ + B as described above. Therefore we can write 

J m 

max D — 2m max ^ — Hj) 2 

Q j=i i=i 


= 2m max || P(AQ + B ) ||p 

r q ' 

— 2m max || (A 0 , B 0 ) 

\_— Ij . 

= V2mmax\\(Ao,B 0 )Q\\ 2 F , 


\l 


where A c and B 0 are centered form of A and B and Q = ^{Q r > —Ij)' G IR^ n+ ‘ / ^ xJ 
and Q'Q = Ij. 

From the above, the form of the minimizing with the cases before is the same 
as that of the shifted version which we are now minimizing the expression of 

min || (A 0 , B 0 )(I n +j - QQ')\\ 2 f = m i n ||A(-A>, B 0 )\\ 2 F . (6.13) 


Thus the results can be summarized as: 
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Correspondence Analysis 

For the shifted version of H = AQ+B, the minimizing problem is equivalent 
to that of 

min \\A(A 0 B 0 )\\ 2 Fl 

or finding the vector X in the orthogonal subspace spanned by Q = , ~h)' 

which satisfies the equation (A 0 , B a )X = 0. 

Total Least Squares 

The total least squares problem AX « B can be stated as follows: find the 
minimum norm vector X = \X\ —2]' which solves exactly the perturbed equation 
(A, B)X = 0, where ( A,B ) = (A, B) — A (A,B). The Frobenius norm of the 
perturbation applied to the appended matrix (A, B), namely ||A(A, has 

to be minimized. Therefore both TLS and CA approaches are equivalent in this 
formulation. 

Assumptions and Observations 

For the TLS and CA equivalence to hold, we have the following assumptions. 

• Both matrix A and B have to be centered as stated above. 

• In Correspondence Analysis, the rank j of the subspace on which the matrix 
A is projected can be chosen freely. However, the rank of the TLS problem 
are usually reduced to the numerical rank of A itself. 
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Chapter 7 

Neural Networks and its 
applications 

7.1 Introduction 

The development of neural networks are motivated by the theory of human neural 
systems. The theory of human neural systems involves the study of the behavior 
of the interconnected processing elements called neurons. Neurons are processing 
or decision making units by which information are passed from one region to an- 
other region. In human nervous system, a vast number of interconnected neurons 
are responsible for processing and relaying information and commands from one 
part of the human body to the other. This complex chain of commands and deci- 
sion making process motivated the development of Artificial Neural Network (NN) 
which mimic the behavior of the human nervous system. The basic elements in 
neural networks called neurons is modeled as that of the human nervous system. 
Signals processed by the neurons travels to adjacent neurons such that a specific 
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command or signal can be passed to the intended destination. In recent years, 
thanks to the availability of high speed computers, much attention has been given 
to the development of Artificial Neural Network with applications ranging from 
pattern recognition, function approximation, image processing, system identifi- 
cation to dynamical system problems. It is illustrated in [20] [12] [21] that a large 
number of problems can be handled by some forms of neural networks without 
prior assumption of the system model. This robust nature makes neural networks 
a good methology for solving a wide range of problems. Before demonstrating the 
usefulness of neural networks in solving our load measurement problem, the fol- 
lowing section will give a brief overview of the theory and development of neural 
networks by presenting some practical examples and applications. 

7.2 An Overview of Neural Networks 

Since the theory of neural networks is still not quite mature, there are still some 
unsolved issues such as the type of architecture of the neural networks, the type 
of training methods and the design procedures. Recently, much works have been 
done on the applications of neural networks. There are two main categories in 
neural networks. The first category is called the dynamic neural networks as the 
system parameters is dependent on time or the network possess ’’memory”. The 
second category is called static neural networks since the system parameters are 
static in time or ’’memoryless”. 
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Figure 7.1: A basic neuron 

7.2.1 Basic Neural Units or Neurons 

Neural Networks consist of basic processing units called neurons. It can be viewed 
as input/output devices where the output is determined by a specific rule or a 
functional relation. The figure (7.1) shows a neuron with multiple inputs x = 
[xi, . . . , x n ]' € IR n and a single output y x . In addition, the neuron is defined by 
its activation function /(•), the bias bj and the weights 
Thus, the output of the i-th neuron can be written as 

n 

Vi = f(Y^ x i w ij +k)> (7-1) 

3 = 1 

where n is the total number of input to the i-th neuron. 

The activation function /(■) is an important element of a neural network for 
the nature of the function introduces nonlinearity into the NN. Such nonlinearites, 
which we will explore later, allows the NN to handle nonlinear problems. The type 
of activation function used is known as the sigmoid function which usually take 
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the form of a step function, tangent, logarithm or sometimes even linear function. 
An examples of such a function is shown in fig (7.2) where the functional relation 
can be expressed as, 


f{x) = 


1 


(7.2) 


1 + e-“ 

Furthermore, the bias 6, and the parameter a determine the shape of the 
sigmoid function which we can select during the design phase of a neural network. 


7.2.2 Basic Neural Networks 

A typical neural network consists of interconnected neurons. Different types of 
architecture are commonly used in NN. Neurons can be fully connected(i.e. each 
neuron is connected to the other neurons via a weigh). However, most popular 
NN requires only neurons to be connected to adjacent layers of neurons which we 
will subsequently discuss in detail. Besides the architecture of the NN, NN can be 
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Input Layer Hidden Layers Output layer 

Figure 7.3: A Multiple Layer Perceptron Neural Network 


classified as the ’’feedforward static NN” and ’’feedback/recurrent /dynamic NN”. 
Feedforward NN allows signals to flow in one direction. One of the important 
types of NN is called Multiple Layer Perceptron (MLP). The second class of NN is 
characterized by the feedback network such as the famous Hopefield Network. [21] 
[32] In chapter 8, NN will be used to solve the load measurement problem. The 
type of NN to be used for this application will be in the form of a static NN or 
MLP. 

7.2.3 Multiple Layer Perceptron (MLP) 

A neural network can be made up of several layers of neurons. In MLP, only 
adjacent layers are connected to each other. Fig. (7.3) shows a fully connected 
MLP Neural Network. 
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A NN can be partitioned into the input layer, the hidden layers, and the 
output layer. In fig (7.3), the number of hidden layers shown can be selected for 
some specific problems. However, the exact number of layers to be chosen for a 
specific problem is still quite arbitrary. It usually requires experimentations or 
trial-and-error methods in determining the number of layers for possible optimal 
solutions. Works in this area shown in [20] proved that two hidden layers NN 
can approximate any types of function if enough number of neurons are being 
used. However, additional layers might reduce the number of required neurons 
and thereby greatly reduces the complexity of the problem. In chapter 9, we will 
discuss this aspect of the problem and suggest a procedure of designing a neural 
network by selecting the number of neurons in a layer using our load estimation 
problem as an example. 

T.2L4 Training Techniques for Neural Networks 

In order for the Neural Network to act as a function approximation, system iden- 
tification or pattern recognition network, we need to train the neural network 
by adjusting the weights and bias iteratively until the desired objective or error 
criterion is achieved. The training procedures can be classified into two types. 
One type is called the supervised learning and the other type is called the un- 
supervised learning. Supervised learning [12] involves the use of target values or 
patterns. Targets and input pairs are presented to the neural network for training 
at each step until an error criterion has been reached. For unsupervised learning, 
the neural network is only presented with inputs without the aid of the targets. 


69 



7.3 Training by Back Propagation (BP) Tech- 
nique 

Due to the relative ease of use and the stability in the solution, Back Propagation 
became one of the most popular training techniques since its introduction [27], 
For enhancing the convergence properties of the BP technique, modifications 
had been made such as the addition of momentum terms in the direction of 
update. BP is essentially based on the optimization theory using first order 
gradient descent methods. The success of the convergence of the BP relies on the 
direction of the search. In addition, the BP attempts to solve the optimization 
problem which is equivalent to the search of the global minimum to an energy 
or error function defined as E. As a result, BP sometimes give solution of local 
minimum and may not always arrive at the global minimum. Other methods 
using a momentum term (or a second order gradient method) such as conjugate 
gradient method as studied in [17] attempts to avoid the solution to settle on a 
local minimum. 

7.3.1 Notations and Definitions 

Here are some of the notations and definitions used in deriving the Back Propa- 
gation method: 

• E - Error or energy function for minimization. 

• L - Last layer of the neural network. 

• Nl~ Number of neurons in the L - th (last) layer. 
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• Ni~ Number of neurons in the Nth layer. 

• m - m-th step of the iterative updating process. 

• ti - Nth element of the training vector. 

• Wij(l ) - The weight connecting j - th neuron’s output in (/ — l)-th layer to i 
-th neurons in Nth layer. 

• a,i(l ) - Output of the Nth neuron in the Nth layer. 

• Ui(l) - Input to the Nth neuron’s activation function in the l layer. 

• (i - Learning rate for the back propagation training. 

• bi(l ) - Bias of the Nth neuron in the Nth layer. 

7.3.2 Derivation 

Using the notations defined previously, the input Ui(l) and output a,(/) to the i 
-th neuron in the l layer can be expressed as 

Ni 

Ui(l) = - i) +6»(0> and 

j - 1 

Oi(l) = /(Uiil)). ( 7 . 3 ) 

In Back Propagation, the weights are updated according to the expression as 

( 7 . 4 ) 

where m is the iteration step and is the change or updates in weights 

for step m. 
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The weight updating process is according to a first order gradient descent 
method such that the updating direction is in the negative gradient of the energy 
or error function E as expressed by 


1 "l , . 


e = ; £ D‘! m ' - 4 m, (i)] 2 . 


(7.5) 


m=l i=l 


where M p is the total number of training patterns and Ni is the number of 
neurons in the last layer(or the number of output elements). 

Using chain rule, the updating step can be expressed as 


a«4”(o = 


= -M: 


dE da\ m) (l) 


da\ m \l) dw K ™’ (/) 


(">)/ 


From equation (7.3) we have, 

da t (m) (/) 
au t (m) (/) 

dutHi) 

dw^\l) 

From equation (7.8) and (7.10), we have 


= /'(Ui m> (0) and 

= oS m) (/-l). 


Au#°(i) = - 1). 

For the last layer, the partied derivative 6j m \L) can be expressed as 


Sj m) (L) = 


dE 


= t ( , m) - a,- m) (L). 


daJ m) (L) 

Again using chain rule the derivative can be expanded as 


(7.6) 

(7.7) 

(7.8) 


(7.9) 

(7.10) 


(7.11) 


(7.12) 


72 



= 


dE 


da 


dE du^il +JQ 

* Ck. (™)/1 . 1\ a^( m )/ 


£J a^ m) (/ + i) aal m; (Q 

Ni+i 


(7.13) 


V 

,(»»»)/ 


= - V aa ; (/ + 1) wi J m) (f + 1) (7.14) 

h da\ m \l + 1) (Z + 1) j 

N| +1 

E + l)/'(«i m> (' + 1))«'5" > (1 + 1),(7.15) 

j=l 

where l — L — 1,. . . ,1. 

Then the weight updating can be written as 

A wfr+'\l) - tog"’ + S m) (0)4 m) (l - 1) 


(7.16) 


Now, 


<5 t (m) (0 

4 m) (L) 


AT < +1 

- E 4"” (I + l)/'(^ ml (i + l))4 m) (i + 1) 

j-i 

(<"> - aS m) (L). 


(7.17) 

(7.18) 


Then we can evaluate, 


Vt-i 


4 m) (t-i) = - e 4“ > (0/'(i4 m, ( i '))”&‘ ) ( t ) 

j=i 

^L -2 

«< m, (L-2) = - e «] m) (i-i)/'(u5 m) (i- i ))4 m ( i - 1 ) ( 7 - 2 °) 


(7.19) 


«< m, (l) = -E«r ) (2)/'(uS m> (2))tvf > (2). (7.21) 

j = 1 

Since can be recursively computed starting from Sj m) (L) = t\ m ^ — 

a- m) (L) at the last layer, the term Back Propagation is used in the sense that the 
error is propagated from the last layer. 
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Update Frequency 


In the last section, the update of the weights are done when each input target 
pair or patterns are presented to the network. This type of update is suitable for 
some on-line problems when real-time training is necessary. However, the update 
can be done when all the patterns have been presented to the network. This type 
of update is called Block Update and is known to be more robust because the 
weights are updated when the training patterns are averaged. Thus for the load 
estimation when real-time training is not necessary, Block Update method will be 
used. The Block Update method can be generalized by the following modification 
to the update equation in equation (7.16) as 

M p 

AwS*'\i ) = «#'(«) + m E T(i - 1 ). ( 7 - 22 ) 

m=l 

where k is the block update step and M v is the total number of training patterns. 
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Chapter 8 

Load Estimation Problem - 
Linear and Non-linear Least 
Squares Approach 


In this chapter, we will discuss the use of nonlinear least square approach in 
solving the load estimation problem. This problem is motivated by the NASA 
F-lll data set. The objective is to use the calibrated data to predict the load 
experienced by the wing structure during the in-flight stage. 

Load Conditions- 

Here we will define the notion of loading conditions in the calibration and in- 
flight stage of the load estimation problem. During calibration, loads are applied 
on different parts of the wing while the corresponding set of n gages’ readings 
are recorded. If we assume that the number of different loading condition is tti, 
we can define the qr-th loading condition as the set of all individual point loads 
as £C (<7 \ where 1 <q < m. 
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Now the individual load applied on the wing structure l = (Zj , . . . ,lj)' € 
IR J(,)x3 can be written as 

li 9) = [«*, SiXi, SiVi], i = J (q) , (8.1) 

where is the relative position of the load on the wing structure and J' q ’ 

is the total number of individual load points for the g-th load condition . 
Single Equivalent Load 

As discussed in chapter 2, the first, second and third elements of the load 
vector li represents the shear, the moment and the torque. In order to simplify 
the set up of the system of linear equations, the notion of equivalent single load 
is used. Thus the equivalent single load e IR lx3 corresponding to the g-th 
load condition CC ^ can be defined as 

L (9) = £Z< (8.2) 

1=1 

= Xl[si,SiZi,s2/i]- ( 8 - 3 ) 

i=l 

For simplicity, we will concentrate on the shear s element of the load vector since 
the use of the other elements will lead to similar procedures and results. Now 
the shear load of CC ^ can be written as 

j(i) 

L (9) = X>. (8.4) 

i=i 

Then the calibration stage can be summarized by the following steps. 
Calibration Stage Procedure 

• Known loads, which are defined as a specific load condition CC ^ = {l\ q ^ \ 1 < 
q < J^}, are placed on different parts of the wing structure as shown in 
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, (q) 

Calibrated Loads of Load Condition CC 



fig. (8.1). Usually, the specific choice of are due to actual loading 
conditions of interest (such as the loads on the leading edge of the wing is 
more pronounced). 


• The magnitude of the known calibrated loads are simultaneously increased (we 
denote as the equivalent load as described in equation (8.4) for 1 < 

j < k where k is the total number of loads of varying magnitude) and the 
corresponding gage readings are recorded (in the Fill Calibration data ex- 
ample, an approximate k = 50 points are recorded). Fig 8.2 shows a typical 
gage vs. load response curve obtained from a specific load condition CC^ q \ 

• Loads are relocated so that a new set of data corresponding to a different 

loading condition can be recorded. 
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• For the Fill data case, there are m = 9 load conditions while each sets 
contains 50 data points. 


As seen in fig. 8.2, high degree of nonlinearities are present in the data. As 
a result, the conventional linear methods are inadequate for accurate estimation. 
However, we will still discuss the linear approach in estimating the load as a 
comparison to the neural network approach. In order to understand the load 
problem, the load estimation procedure can be summarized in fig (8.3). 
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Figure 8.3: The Load Estimation Problem 



8.1 Linear Approach 

The conventional linear approach attempts to find a relationship between the 
gage and load data (Mb ~ L) in a linear least square sense. It solves for the 
minimizer b of the system of linear equations such that 


mm\\Mb - L\\ = \\Mb - L\\. 

b 


(8.5) 


Since the gage matrix M G IR mfcxn and the load vector L G the system of 

equations can be written as 

Mb « L 
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m is the number of loading conditions, 


( 8 . 6 ) 


(8.7) 
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k is the number of points per loading condition, and 
n is the number of gages. 

The solution b to the linear system can be expressed as 

b = M+L, (8.8) 

where M + is the Moore-Penrose pseudo inverse of M. 

The corresponding load estimate L can then be computed as 

L = fh'M + L > (8.9) 

where fh is a n x 1 vector from the in-flight gage reading. 

8.2 Non Linear Least Squares Approach 

Since nonlinearities are present in the data set, the load problem of solving the 
nonlinear system(/(M) « L) can be written as 

min \\f{M) — L\\, (8.10) 

0 

where f(M, b ) is a nonlinear function of M and 6. 

8.2.1 Neural Network Approach 

In solving the problem min*, ||/(M) — L\\ where /(•) is a nonlinear function, the 
neural network approach can be used as an approximator of the function /(•). 

First, the structure of the NN has to be chosen for this problem. Since there 
is no definite approach in choosing the number of layers or the number of neurons 
for the network, we need to use a trail-and-error approach to determine a possible 
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Targets 


T 



Figure 8.4: The Back Propagation Training Process 

optimal selection. In chapter 9, we will discuss some aspects of determining the 
number of neurons in a hidden layer. In our simulations, two hidden layers are 
used for it is found to be adequate in most cases. The Back Propagation (BP) 
training procedure, which can be illustrated in fig. (8.4), is used in the learning 
process. 

Once the Network is properly trained (i.e. when an error criterion has been 
reajched), the NN can be used to estimate the wing load during the in-flight stage. 
The prediction process is fast enough for real-time application since all the time 
consuming computation or training can be done in the ground calibration phase. 
Fig (8.5) shows that a load estimate can be obtained by passing the in-flight gage 
data to the trained Neural Network. 
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Figure 8.5: The Load Prediction by a Trained Neural Network 

8.3 Simulation Results and Observations 


In order to evaluate the performance of the trained Neural Network, the NN has 
to be tested using existing calibration data. The sigmoid function used for this 
example is defined in (7.2), where we set a — 1. Although we have m = 9 sets 
of loading conditions, we only used m = 5 sets for calibration and training. The 
remaining load conditions are used to validate and quantify the performance of 
the network. After the training is completed, the NN is tested with three different 
sets of loading conditions. In this simulation, the performance of the NN and the 
linear methods are compared by computing the relative error which is denoted 


by 



( 8 . 11 ) 


and similarly the relative error of the linear approach can be denoted as 


£ls 


Ltrue ^LS, estimate 

Ltrue 


II. 


( 8 , 12 ) 


where 

L true is the true calibrated load used as a verification, 

Lnn, estimate is the estimated load using Neural Network approach, and 
Lls, estimate is the estimated load using Neural Network approach. 
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Three performance analysis plots in fig. 8.6, 8.7 and fig. 8.8 show e/v/v and e^s 
vs. load test samples. 
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Relative Error 



Test Load Samples - Load Condition 1 

Figure 8.6: Relative Error vs. Test Load Samples for NN and LS approach(Load 
Condition 1) 


Q'J-. 



Relative Error 



Test Load Samples - Load Condition 2 

Figure 8.7: Relative Error vs. Test Load Samples for NN and LS approach (Load 
Condition 2) 
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Relative Error 


0.5 


* Linear LS method 
o NN method: 



Test Load Samples - Load Condition 3 

Figure 8.8: Relative Error vs. Test Load Samples for N and LS approach(Load 
Condition 3) 


As observed in fig. 8. 6,8.7, and fig. 8.8, the simulation results from the Neural 
Network approach* have better performance in most sample cases. In addition 
large errors occurred in using linear LS method for some specific data samples, 
whereas the errors’ from the NN Approach are all relatively small. In general, 
the relative error ofi NN are more or less confined in the range of 20 ~ 30% or 
below. In addition?,, for the linear LS method, large error occured in some specific 
examples. This seems to be a verification of the general superiority of the NN 
to the LS method when nonlinearlities are present in the system. Moreover, NN 
approach performance will not deteriorate even when the data is actually linear 
for NN can handle both linear and nonlinear systems. 
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8.4 Distributed Load Estimation Problem 


In the past section 8.2.1, the solution for the single load estimate are presented. 
Single load estimate is the equivalent load estimate which consists of a single 
load point L = (S, Sx, Sy), where S denote the equivalent shear load and ( Xi,yi ) 
denote the relative location of the equivalent load on the wing surface. It provides 
only the location and the magnitude of the equivalent load(or equivalently the 
shear, moment and torque components of the Load L). In real life situations, 
loads are usually distributed on different locations of the wing structure. Thus 
the knowledge of the single equivalent point load is not adequate in providing 
information as how the load is distributed on the wing. This lack of information 
might becomes undesirable in dangerous situations such as when a large load is 
concentrated on a specific region of the wing structure. As a result, the estimate 
of the distributed load is highly relevant and useful since it might provide an 
’’early warning signal” for the presence of structure overloading. 

Fig. 8.9 shows the difference between the single load and the distributed load 
representation. It also shows how the wing structure can be partitioned into 
regions of interest. 

Distributed load estimate not only provides the estimate of the single equiv- 
alent load, but also provides the positions and magnitudes of all the forces in 
specific regions of the wing. 

Distributed Load Problem The distributed load problem selects the right pattern 
q = q* or CC^ = < j < J^} from the set (CC^\CC ^ 2 \. . . 

which represents the pattern set of all possible flight conditions. 

If we can solve this general problem, we will be able to identify the magnitude 
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Figure 8.9: Single Equivalent Load and Distributed Load on a Wing Surface 
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of the specific force acting on a specific small region of the wing. As a result, the 
solution might provide valid information such as early warning signals for local 
overloading of the wing structure. However, such a task seems to be too ambi- 
tious. Alternatively, a sub-optimal approach can be adopted in order to retrieve 
as much information as possible for the in-flight data analysis. Since we only have 
a limited number of load patterns obtained in the calibration stage, the best we 
can do is to utilize those sets of data to achieve our goal. If each patterns in 
the calibration stage are valid in-flight characteristic loading conditions, we can 
use them as the training patterns. Once the neural network has been properly 
trained, we could(with some degree of accuracy) be able to identify which previ- 
ously trained load pattern (or loading condition) the wing is experiencing during 
the in-flight session. Such a procedure seems useful since if we train the network 
with some specific conditions that represent a dangerous overloading situations, 
’’early warning signals” might be available once this pattern is encountered during 
flight. Before apply the neural networks approach to the distributed load prob- 
lem, the next section will discuss some basics of pattern recognition by neural 
networks. 


8.5 Pattern Recognition by Neural Network 

Pattern recognition by neural networks can be described as using nonlinear clas- 
sifier to separate classes or patterns by drawing partitions between classes in the 
pattern space. Fig (8.10) shows a typical pattern space and how nonlinear sepa- 
rable cases can be partitioned into decision regions using nonlinear classifier such 
as NN. 
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The training input matrix and target matrix for this neural network are the 
gage matrix M and the target matrix T where 
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and 


(8.13) 


(8.14) 


The neural network approach in solving the problem of identifying the pat- 
terns of loading can be summarized as follows: 


• The structure of the neural network such as the number of layers and the 
number of neurons are chosen by in some(ad-hoc) maimer (chapter 9 will 
discuss the issue of the selection of the number of neurons). 


• M' and T are presented to the neural network and trained by back propa- 
gation (using block update frequency as described in the last chapter) until 
the maximum number of steps has exceeded or the the error criteria e/ is 
is met as given by 

E — \\T — TH < c/, (8.15) 

where T is the output of the neural network. 

• A noisy version of the gage matrix which we define as M and the target 
T are presented to the Neural Network. Then the second phase of training 
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continues until the error criteria is met or the maximum number of steps 
has been exceeded. The noisy version of the gage matrix can be expressed 


M> = 


(8.16) 

itf'i = 

[M' + e j], and 

(8.17) 

M' 2 = 

[M' + e'], 

(8.18) 


where e\ and e 2 are mkx n normal random noise matrix with two different 
amount of noise variance. 

• The Neural Network is again trained with the noiseless matrix M' and T 
to guarantee accurate prediction of the patterns for noiseless cases. 

• Once the NN has been trained, the in-flight gage data fh will be passed into 
the NN. A ’’one” in the z-th position of the output target f indicates the 
presence of the z-th load pattern. 


8.6 Results and Observation 

Results shown in fig. (8.11) and. (8.12) shows reasonably good recognition per- 
formance even for noisy situations. For 10% noise level, a recognition error of 
less than 1% is observed- 
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2 hidden layer NN for 2 classes recognition with 20, 2 neurons in layer 1, 2 
Figure 8.11: Neural Network for Load Pattern Recognition 
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Noise Level 


2 hidden layer NN for 3 classes recognition and 24, 4 neurons in layers 1, 2 


Percentage of Recognition Errors 



Noise Level 

2 hidden layer NN for 4 classes recognition 32, 8 neurons in layer si, 2 
Figure 8.12: Neural Network for Load Pattern Recognition 
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Chapter 9 

Performance Analysis and 
Design Criteria of Neural 
Networks 


In this chapter, we will discuss the issue of design criteria and the performance 
aspects of neural networks. There are a number of parameters determining the 
performance (accuracy) of the neural networks in estimating the targets(loads). 
The parameters can be listed as follows. 

• e = error criteria for stopping the training process. 

• n = number of iteration required for convergence. 

• p = number of steps or time needed per iteration n. 

• a e = the noise variance for the model. 

• a = the degree of nonlinearities in the sigmoid function. 
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• rid = total number of nodes in the neural network. 


• rii = number of hidden layers in the neural network. 

The two basic important terms which quantify the performance of the neural 
network can be defined as the prediction accuracy p y and the time t required for 
training. 

• t = n* p = total time needed for the training process. 

• pj = the load prediction accuracy. 

The total time for training can be related as a function of the basic parameters 
as 

t = f(e,(T e ,a,n d ,ni). (9.1) 

Similarly, the performance of the neural network can also be related as 

p f - g(e,<r e ,a,n d ,ni). (9.2) 

Fig. 9.1 shows the simulation results of the F-lll dataset(the gage and load 
data are listed in Appendix B) in which the number of neurons in the first hidden 
layer is varied. The number of neurons in the second layer are fixed at 2. Since 
excessive neurons often causes overfitting of the model and insufficient neurons 
causes underfitting of the model, we expect that the errors will be high in these 
cases. Indeed from Fig. 9.1, we observed that the relative errors are both high in 
these cases. Thus we can conclude from this example, the optimum choice seems 
to lie in the vicinity of 14 neurons in the first hidden layer. 
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Average Relative Error 




9.0. 1 Generalization and Training Accuracy 

The ability of the NN to estimate the targets with untrained data is sometimes 
referred to the term generalization performance [17]. The error criteria e can also 
be interpreted as the training accuracy of the NN. Thus the performance of the 
NN depends more on the generalization performance rather than the training 
accuracy of the network. The issue of generalization and training accuracy can 
be viewed as an analogy to the issue of noise sensitivity and model accuracy as 
discussed in chapter 5. Although a large number of neurons can increase the 
training accuracy, the problem of overfitting might cause the network to have 
a poor generalization performance. Thus we need to select the right number of 
neurons such that the effects of overfitting can be avoided or reduced. In other 
words, we need to increase the number of neurons only if it is necessary. 

In order to quantify the generalization performance of the network, an ap- 
proach called ’’Cross Validation”, in which calibration data are separated into 
two parts, can be used. One part of the data are used for training and the other 
part are used for validation. Thus from the validation results, one can quantify 
the generalization performance of the network. 

9.0. 2 SVD, CA and Collinearity Method for Reducing 

Neurons in a Layer 

Recently, Xue in [29] proposed the use of SVD in determining the number of 
neurons in a hidden layer. By determining the effective rank of the weight matrix, 
one can effectively determine the number o£ neurons needed in a hidden layer. 

Let the output of the z'-th hidden layer be A £ IR m , the input of the i-th layer 
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be P G IR n , where m is total number of neurons in the i-th hidden layer, and n 
is the total number of inputs from the z-th hidden layer, and W € IR mxn is the 
weight matrix of the z-th hidden layer. 

Thus the output of the z-th layer can be written as, 

A = WP, (9-3) 

where each row of W represent the weights of a neuron. If two rows of W are 
nearly identical, the output vector A will have two identical elements. Thus if 
two neurons give the same output for every input vector A , reducing one neuron 
from these two will give the same results. Based on this observation, if two rows 
of W are nearly linearly dependent, one can eliminate the extra neuron without 
affecting the performance of the network. 

However in practice, the determination of the effective rank of a matrix is 
usually a sensitive issue. Xue in [29] proposed the use of an error criteria in 
determining the effective rank of the weight matrix W . However, the use of this 
method is not straightforward and clear. By using correspondence analysis(CA) 
and collinearity index(CO), we can provide an alternate and simple method to 
determine the rank of W and thus reduce the number of neurons in a layer. 

Correspondence analysis and collinarity method are applied to our simulation 
example as shown in fig. (9.1). From fig. (9.2) and fig. (9.3) we can group the 
nearby point to form one single point. Using CA, the effective rank of the weight 
matrix IF is 15 while using CO gives an effective rank of 13. Thus this simple 
example seems to agree with our simulation shown in fig. (9.1) that the vicinity 
of 16 neurons is in the region of the optimum choice of neurons. Although the use 
of the CA and CO is not straightforward as the grouping of points requires the 
right criteria of closeness. Different criteria of closeness between points sometimes 
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first dominant sing, vector 

Figure 9.2: Correspondence Analysis of the weight matrix W 


render different results. However we provided a practical and simple approach 
of combining SVD, CA, and CO in determining the possible optimal number of 
neurons in a layer. 
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Chapter 10 

Conclusions and Future Work 


This thesis provided a study of parameter estimation using linear and nonlinear 
LS approach. We have discussed the issue of model accuracy /noise sensitivity 
and derived a bound for the change in the model accuracy (represented by its 
residual) as a function of the noise. The problem of the trade-off between model 
accuracy and noise sensitive is a major issue. This problem is also analogous 
to the issue of generalization and training accuracy of the neural network. We 
have also studied and derived the equivalency of TLS and CA. These two seem- 
ingly different topics in the field of signal processing and statistical analysis are 
shown to be equivalent with proper scaling and centering of the data matrix. 
In practice, the data matrix are sometimes nonlinear in nature and the use of 
linear LS methods is insufficient for accurate parameter estimations. Motivated 
by the nonlinearities in the NASA F-lll data matrix, we studied the possible 
use of nonlinear LS methods. Conventional nonlinear least square methods and 
spline approximation require the exact modeling of the functional relationship. 
The use of neural network provides a practical and robust approach in dealing 
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with highly nonlinear data. Using the NASA data as the training data, we sim- 
ulated the estimation procedure by using neural network approaches. Although 
the computational effort in using neural networks is very high, the results we ob- 
tained are both useful and encouraging. Compared to the conventional linear LS 
approach, the neural network approach yields significantly better performance in 
load estimation. We also showed the ability of neural network in estimating the 
load conditions. The recognition of the trained load conditions during flight pro- 
vided a possible method to provide ” early warning signal” for the wing structure 
once a dangerous overloading condition is encountered. 

We also studied the aspects of designing neural networks. The selection of 
the number of neurons in a hidden layer is crucial in designing a neural network 
with good generalization and estimation accuracy. Underfitting and overfitting 
yields poor performing network if the number of neurons is not correctly used. 
Work done by Xue in [29] proposed the use of effective rank of the weight matrix 
determined by SVD as a way to reduce the number of neurons in a layer. We 
found that Correspondence Analysis and Collinearity Index methods appeared 
to be equally effective in determining the effective rank of the weight matrix and 
thereby effectively reduces the number of neurons in a hidden layer. 

We have studied both the linear and nonlinear methods for solving LS prob- 
lem. As we expected, the linear methods are more analytically tractable while 
nonlinear methods need to be based on heuristic techniques. Thus more works is 
needed to fully understand the nonlinear methods including the neural network 
approaches. 

Some of the unsolved problem and possible extension to our work include: 

• Other types of NN other than MLP can be used for the Load Problem. 
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• Analytical analysis and design aspects of NN need to be considered. 

• The use of a mix of NN types to recognize complex decision regions such 
as the ’’exclusive or(XOR)” operation (disconnected pattern space), which is 
not possible for MLP type NN. 

• The use of basis function in defining the pressure distribution can be used 
instead of point loads. 

• Computer simulated structural analysis can simulate gage-load relation but 
lack the ability to model the actual structure itself. 

• A combination of computer analysis and structural testing(calibration) can 
be used to design a better and more accurate flight load estimations system. 
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Appendix A 

Computer Codes for Neural 
Network Load Estimation 


This appendix lists the matlab [6] computer codes for the single load estimation 
procedure. The procedures are summarized below: 

1. datreduct.m Reduces the number of points in a load condition; 

2. datred.m Select the load conditions to be reduced; 

3. datsetup.m Set up the gage/load data for training; 

4. trainNW.m Ttain the NN with noise; 

5. trainNN5.m Train the NN without noise; 

6. testset.m Test the NN performance. 
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datreduct .m 

res = 1; 

while (res==l) f 

clg; 

hold off; 

Lef f =b; 
mm=a ; 

ss=size (a) ;ss=ss (1) /ss^ss/S-S/ss^roundtss) ; 

Leff=Leff (l:ss) ; 

mm=mm(l : ss, : ) ; 

plot (Leff , mm, 1 * 1 ) ; 

hold on; 

plot (Lef f , mm) ; 

s=size (Lef f ) ; s=s (1) ; 

step= (max (Lef f ) -min (Lef f ) ) /s (1) ; 

%npoint=input ( ’please input number of points 1 ); 

disp ('the maximum number of points you can input is ' ); 
disp (s) ; 

disp{ ’please input the Leff value to be picked for your reduced 
data ’ ) ; 

disp(’e.g. input [-2000 -4000 -5000 -8000 -11000] for 5 points, 
but more than 1 point ’ ) ; 

disp (’or the name of the vector, e.g. v ’); 

Lin=input ( ’please enter now ’); 
np=size (Lin 1 ) ; 
npoint=np (1) ; 

mmr=zeros (npoint , 12) ; Lef f r=zeros (npoint, 1) ; 
for k“l: npoint, 
temp « 1; 
for i=l:ss, 

tempi = Leff(i) - Lin(k); 
if sign (tempi) — -sign (temp), 
mmr (k, : ) =mm(i, : ) ; 

Leffr (k)=Leff (i) ; 
temp Leff(i) - Lin(k); 

% disp (temp); 

else 

temp = Leff(i) - Lin(k); 

end; 

end; 

end; 

plot (Leff r, mmr) ; 
plot (Lef fr, mmr, ’o’); 
hold off; 
gagemat=mmr ; 
gagemat=-gagemat * ; 

% gagematrix normalization 
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mm=max (max (gagemat) ) ; 


gagemat =gagemat /mm; 
targets=Lef f (i) ; 

% targets normalization 
targets --targets ' ; 
tt==max (targets 1 ) ; 
targets=targets/tt ; 

disp( f enter 1 if you want to fit better ’); 
res~input ( ’enter 0 if you are satified »’); 

end. 
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datred.m 

disp( , load condition 1.3 LC 1 ') 

disp('load condition 11.3 LC 2 f ) 

disp('load condition 13.3 LC 3 ’) 

disp(’load condition 17.2 LC 4 ') 

disp(’load condition 21.4 LC 5 1 ) 

disp(’load condition 24.8 — — LC 6 1 ) 

disp(’load condition 28.4 LC 7 ’) 

disp(*load condition 9.3 LC 8 1 ) 

disp(’load condition 91011 LC 9 1 ) 


t <1, : ) «’dat01 1 ; 
t <2, : ) * 'dat 02 1 ; 
t (3, : ) = 'dat 03 ' ; 
t (4, : ) = 1 dat 0 4'; 
t <5, : ) = ' dat 05'; 
t <6, :)«'dat06' ; 
t (7, : ) = ' dat 07 ' ; 
t (8, : ) = 1 dat 08 ' ; 
t (9, : ) = 1 dat 09 1 ; 


tr (1, : ) « f rdatOl f 
tr (2, : ) = ' rdat02 ' 
tr (3 f : ) -’ rdat03 1 
tr ( 4 f : ) «= 1 rdat 04 1 
tr (5, : ) = ? rdat05 1 
tr (6, : ) = 1 rdat 06 1 
tr (7, : ) = f rdat07 1 
tr ( 8 r : ) = 1 rdat 08 1 
tr ( 9 , : ) = 1 rdat 0 9 f 


disp ( ’enter the load conditions to be used for training *) 
dispt’For example [12345] for selecting LC1 to LC5 r ); 
kk-input ( ’enter 1 ) ; 


s^size (kk) ; s=s (2) ; 


for i*l:s f 
j-kk (i) ; 
eval ( t { j , : ) ) ; 
date reduct ; 

c=setstr { [ ' rdat 1 , int2str ( j) ] ) ; 
if j“l r 

save rdat 01 mmr /ascii /double; 
save rdat01_L Leffr /ascii /double; 
elseif j==2,; 

save rdat 02 mmr /ascii /double; 
save rdat02_L Leffr /ascii /double 
elseif j“3, 

save rdat 03 mmr /ascii /double; 
save rdat03_L Leffr /ascii /double 
elseif j“4, 

save rdat 04 mmr /ascii /double; 
save rdat04_L Leffr /ascii /double 
elseif j ——5 r 

save rdat 05 mmr /ascii /double; 
save rdat05_L Leffr /ascii /double 
elseif j==6, 

save rdat 06 mmr /ascii /double; 
save rdat06 L Leffr /ascii /double 
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elseif j“7, 

save rdat07 mmr /ascii /double; 
save rdat07_L Leffr /ascii /double 
elseif j“8, 

save rdat08 ininr /ascii /double; 
save rdat08_L Leffr /ascii /double 
elseif j==9, 

save rdat09 mmr /ascii /double; 
save rdat09_L Leffr /ascii /double 
end; 

end; 

end. 
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dat setup . m 

disp{' Enter the number Load conditions to be included for 
training 1 ) ; 

disp('For example, enter [12345] LC1 to LC 5 is used for 

training' ) ; 

temp3=input ( ' ent er» ' ) ; 

temp4=size (temp3) ; 

temp4«temp4 (2) ; 

gmmr- [ ] ; 

gLef f r= [] ; 

for i * l:temp4, 

ttt^tempS (i) ; 

if ttt“l, 

load rdatOl; 

load rdat01_L; 

mmr= rda 1 01; Le f f r *= rda 1 0 1_L ; 

elseif ttt“2,; 

load rdat02; 

load rdat02_L; 

mmr=rdat02 ; Lef f r=rdat02_L; 

elseif ttt"3, 

load rdat03; 

load rdat03_L; 

mmr = r da 1 0 3 ; Le f f r = rda 1 0 3_L ; 

elseif ttt==4, 

load rdat04; 

load rdat04_L; 

mmr* rda 1 04; Le f f r = r da 1 0 4_L ; 

elseif ttt~5, 

load rdatOS; 

load rdat05_L; 

mmr* r da 1 05; Le f f r = r da 1 0 5_L ; 

elseif ttt==6, 

load rdat06; 

load rdat06_L; 

mmr=rdat06; Lef f r=rdat06_L; 

elseif ttt***7, 

load rdat07; 

load rdat07_L; 

mmr=rdat07 ; Lef f r=rdat 07_L; 

elseif ttt==8, 

load rdat08; 

load rdat08_L; 

mmr = r da 1 0 8 ; Le f f r = r da 1 0 8_L ; 

elseif ttt«==9, 

load rdat09; 

load rdat09_L; 

mmr = rdat09; Leffr=rdatO 9_L ; 

end ; 

gmmr = [ gmmr ' , mmr ' ] 1 ; 
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gLef f r= [gLef f r 1 ,Leffr T ] 1 ; 
end; 
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trainNW . m 


% Train pattern recognition problem 

% load conditions are represented as patterns 

% 

% 

% TRAIN WITH NOISE 
% 

% This file trains a two hidden layer log- sigmoid/ log- sigmoid 

% network to classify loading conditions patterns 

% 

% 

% The network is first trained on noise free loads. 

% It is then trained on loads with noise. 

% the network is then again trained on noise free loads. 

% The result is a network which can properly estimate 

% noise free loads and does a good job of estimating 

% loads with noise. 

% LOAD PROBLEM 

% 

% INITIALIZE NETWORK ARCHITECTURE 

[R,Q] = size (gagemat) ; 

51 =20; 

52 =2; 

disp ('the number of 1 st hidden layer neurons is '); 
disp (SI) ; 

res=input ( 1 enter y if you want to change it/ otherwise push enter 
» 1 , 1 s ' ) ; 
if res=='y', 

Sl=input ( 'enter no of 1st layer neurons »'); 
end; 

disp ('the number of 2 nd hidden layer neurons is ') ; 
disp (S2) ; 

res=input ( 'enter y if you want to change it, otherwise push enter 
» ' , ' s ' ) ; 
if res=='y'/ 

S2=input ( 'enter no of 2nd layer neurons »'); 
end; 

[S3, Q] = size (targets) ; 

[W1,B1] = nwlog(Sl,R); 

W2 = rands(S2,Sl) *0.01; 

B2 = rands (S2, 1) *0 .01; 

W3=rands (S3,S2) *.01; 

B3=rands (S3, 1) * . 01; 

% TRAIN THE FUNCTION 
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% PROBLEM 

% Get gage data to be trained 
P *= gmmr 1 ; 

T = gLef f r * ; 

% Normalize the data 
nmm=max (max (abs (P) ) ) ; 
ntt~max (abs (T) ) ; 

P=-P . /nmm; 

T— T/ntt; 

% TRAINING PARAMETERS 
disp_freq * 200; 
max__epoch * 16000; 
err_goal = 0.01; 
lr * 0.05; 
lr_inc * 1.05; 
lr_dec = 0.7; 
momentum - 0.95; 
err_ratio *= 1.04; 

% TRAINING WITHOUT NOISE 

TP * [disp_freq max_epoch err_goal lr lr_inc lr_dec momentum 
err_ratio] ; 

[W1,B1,W2, B2,W3,B3,xx]=trainbpx(Wl,Bl, ’logsig' ,W2,B2, ... 

T logsig’ , W3, B3, 1 logsig’ , P, T, TP) ; 

% SAVE NETWORK TRAINED WITHOUT NOISE 

save prl_wl W1 /ascii /double 
save prl_bl B1 /ascii /double 
save prl_w2 W2 /ascii /double 
save prl_b2 B2 /ascii /double 
save prl_w3 W3 /ascii /double 
save prl_b3 B3 /ascii /double 
% TRAINING PARAMETERS 
max_epoch ■= 3000; 
err_goal = 0.06; 
rand ( ’uniform’ ) ; 

TP «= [disp__freq max_epoch err_goal lr lr_inc lr__dec momentum 
err_ratio] ; 

% TRAINING WITH NOISE 
for pass - 1:5 

fprintf ( ’Pass = % . Of \n 1 , pass) ; 

[ddl, dd2 ] -size (gagemat ) ; 

%avegagval s =norm (gagemat) /ddl / dd2 ; 

P * [gagemat, gagemat, . . . 

(gagemat + gagemat . *rand(R, Q) *0 . 1) , ... 

(gagemat + gagemat . *rand (R, Q) *0 . 2) ] ; 

T = [targets targets targets targets] ; 

%size (T) ; 

[Wl,Bl,W2,B2,W3,B3,xx]=trainbpx(Wl,Bl, ’ logsig ’, W2 , B2, . . . 
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'logsig* ,W3,B3, * logsig’ , P, T, TP) ; 
end 

rand ( 'uniform’ ) ; 

% PROBLEM 

P ■ gagemat ; 

T * targets; 

% TRAINING PARAMETERS 

max_epoch » 16000; 
err_goal = 0.05; 

% TRAINING WITHOUT NOISE AGAIN 

TP = [disp_freq max_epoch err_goal lr lr_inc lr_dec momentum 
err_ratio] ; 

[W1,B1,W2, B2,W3, B3, xx] ^trainbpx (Wl, Bl, 1 logsig 1 , W2, B2, . . . 
'logsig' , W3, B3, 1 logsig' ,P,T, TP) ; 

% SAVE NETWORK TRAINED WITH NOISE 

save pr2_wl W1 /ascii /double 
save pr2_bl B1 /ascii /double 
save pr2_w2 W2 /ascii /double 
save pr2_b2 B2 /ascii /double 
save pr2_w3 W3 /ascii /double 
save pr2_b3 B3 /ascii /double 
% SUMMARIZE RESULTS 

A - logsig (w3*logsig (W2*logsig (W1*P, Bl) , B2) ,B3) ; 

SSE * sumsqr(A-T); 

fprintf { ’Final sum squared error without noise: %g.\n',SSE); 
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trainNNS . m 


% Train Neural Networks problem 

% load conditions are represented as patterns 

% Using gage matrix of datset 

% This file trains a two hidden layers log-sigmoid/ log-sigmoid 

% network to classify loading conditions patterns 

% 

% 

% The network is first trained on noise free loads. 

% The result is a network which can properly classify 

% noise free loading conditions and does a good job of 

classifying 

% loads with noise . 

% LOAD PROBLEM 

% 

% INITIALIZE NETWORK ARCHITECTURE 

[R,Q] - size (gagemat) ; 

51 -20; 

52 =2; 

disp('the number of 1 st hidden layer neurons is f ); 
disp(Sl); 

res-input ( 1 enter y if you want to change it, otherwise push enter 
» 1 , ' s 1 ) ; 
if res== 1 y ' / 

Sl=input ( ’enter no of 1st layer neurons »’); 
end; 

dispCthe number of 2 nd hidden layer neurons is '); 
disp(S2); 

res=input (’ enter y if you want to change it, otherwise push enter 
» ’ , ’ s ' ) ; 
if res— ’y’, 

S2-input ( ’enter no of 2nd layer neurons »’); 
end; 

[S3, Q] - size (targets) ; 

[W1,B1] = nwlog (SI, R) ; 

W2 = rands (S2, SI) *0.01; 

B2 = rands (S2, 1) *0 . 01; 

W3=rands (S3, S2) *.01; 

B3= rands (S3, 1) * . 01; 

% TRAIN THE FUNCTION 

% PROBLEM 

% Get gage data to be trained 
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P = gmmr ' ; 

T « gLef f r 1 ; 

% Normalize the data 
nmm=max (max (abs (P) ) ) ; 
ntt^ax (abs (T) ) ; 

P*=-P . /nmm; 

T— T/ntt; 

% TRAINING PARAMETERS 
disp_freq = 10; 
max_epoch * 16000; 
err_goal *= 0.005; 
lr - 0.05; 
lr_inc * 1.05; 
lr_dec - 0.7; 
momentum = 0.95; 
err_ratio * 1.04; 

% TRAINING WITHOUT NOISE 

TP * [disp_freq max_epoch err_goal lr lr_inc lr_dec momentum. 
err_ratio] ; 

[Wl, Bl, W2, B2, W3,B3,xx]=trainbpx(Wl, Bl, 1 logsig' , W2, B2 , . . . 

1 logsig' , W3, B3, 1 logsig 1 , P, T, TP) ; 

% SAVE NETWORK TRAINED WITHOUT NOISE 

save prl_wl Wl /ascii /double 
save prl_bl Bl /ascii /double 
save prl_w2 W2 /ascii /double 
save prl_b2 B2 /ascii /double 
save prl__w3 W3 /ascii /double 
save prl_b3 B3 /ascii /double 
% TRAINING PARAMETERS 

% SUMMARIZE RESULTS 

A * logsig (W3* logsig (W2* logsig <W1*P, Bl) ,B2) , B3) ; 

SSE = sumsqr(A-T); 

fprintf ( 'Final sum squared error without noise: %g.\n',SSE>; 
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testset.m 

% This is to test the NN and Normal equation approach 
% For load condition 1 to 9 the relative error is displayed 
% The relative error is stored in errn= [error-Normal 
Equation, error-NN, error-NN-t rained with noise] 

% error-NE relative error using normal equation approach 

% error-NN relative error using NN and trained without noise 

% error-NN-WN relative error using NN and trained with noise 


disp ('enter the Load condtion to be tested 1 ); 
in=input ( ' enter» ' ) ; 

for k=l : 9 
if k-=in, 

eval ( [ 'datO 1 , int2str (in) ] ) ; 

end; 

end; 

disp( ? enter 1 if the NN is trained with noise'); 
disp ('enter 0 if the NN is trained without noise'); 
temp9=input ( ' enter» 1 ) ; 


% Test with no noise trained network 

load prl__wl; Wl=prl_wl; 

load prl_bl; Bl*prl_bl; 

load prl_w2; W2=prl__w2; 

load prl_b2; B2=prl_b2; 

load prl_w3; W3=prl_w3; 

load prl_b3; B3=prl_b3; 

if temp 9 == 1, 
load pr2_wl; 
load pr2_bl; 
load pr2_w2; 
load pr2_b2; 
load pr2__w3; 
load pr2_b3; 
end; 

% Linear solutions 
AA=P ' ; 

x=pinv (AA) *T 1 ; 
res=norm(AA*x-T 1 ) ; 
res=res*res; 
szl=size (b) ; 
szl=szl (1) ; 
szl=round (szl (1) /2) ; 
a=a (1 : szl, : ) ; 
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b=b (1 : szl, : ) ; 

%a=a (1 :50, : ) ; 
%b-b(l:50, :) ; 


Pt*=-a 1 /mrni; 
targetss--b' /ntt ; 


%P=n2gf 11 ; 
%targetss=n2Lll; 


% A is the output of the NN network — or estimated load using NN 
% P is the test gage values 

% T is the targets or true load for testing 


A = logsig (W3*logsig (W2*logsig (Wl*Pt, Bl) , B2) , B3) ; 
if temp 9 —1, 

A2 * 

logsig (pr2_w3*logsig (pr2_w2*logsig (pr2_wl*Pt, pr2_bl) ,pr2_b2) , pr2__b 

3); 

end; 

error=A-targetss ; 

YY=Pt 1 *x; 

errorN=YY T -t arget ss ; 

errorNn=errorN . /t arget ss; 
errorn=error. /target ss; 
errn= [errorNn* ,errorn’ ] ; 

if temp9“l, 

er ror2=A2-t arget ss ; 

errorn2=error2 . /target ss; 

errn= [errorNn 1 , errorn * , errorn2 * ] ; 

end; 

%plot (abs (error) , ' * * ) ; 

plot (abs (errn) , 1 * ' ) ; 

grid; 

end; 
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Appendix B 

Calibration Data from NASA 
F-lll Load Measurement 


The reduced load vector L G IR 50 and the corresponding gage matrix M G 
IR 50x12 are listed below. 
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- 3,571 v 
— 5 , 561 
- 4,710 
- 7,010 
- 2,833 
- 3,544 
- 2,087 
- 2,653 
- 3,554 
- 4,571 
- 1,727 
- 1,755 
- 2,073 
- 4,191 
- 5,314 
- 7,010 
- 0,795 
- 11,669 
- 1 , 704 
- 1 , 710 
- 4 , 110 
- 5 , 210 
- 5,952 
- 8, 565 
- 11,190 
- 13,116 
-671 
- 1,397 
- 2,293 
- 3 , 101 
- 3,887 
- 4,660 
- 6 , 129 
- 7,045 
-301 
-384 
-779 
- 1,718 
- 2 , 465 
- 3,216 
- 4 , 721 
- 5,655 
-739 
-737 
- 1 , 173 
- 2,870 
- 4,240 
- 5,587 
- 8,317 
- 10 , 025 / 


— 0.23 

- 0.26 

- 0.12 

- 0.12 

- 0.45 

- 0.60 

- 0.45 

- 0.63 

- 0.69 

- 0.75 

- 0.75 

— 0 . 88 \ 

- 0.52 

- 0.46 

- 0.24 

- 0.23 

- 0.96 

- 1.25 

- 0.91 

- 1.22 

- 1.45 

- 1.54 

- 1.46 

- 1.71 ' 

- 0.34 

- 0.26 

- 0.16 

- 0.15 

- 0.66 

- 0.85 

- 0.59 

- 0.78 

- 0.99 

- 1.03 

- 0.96 

- 1.13 

- 0.82 

- 0.78 

- 0.39 

- 0.39 

- 1.56 

- 2.02 

- 1.47 

- 1.99 

- 2.33 

- 2.49 

- 2.36 

- 2.75 

- 0.32 

- 0.18 

- 0.15 

- 0.11 

- 0.67 

- 0.91 

- 0.53 

- 0.76 

- 1,11 

- 1.13 

- 0.98 

- 1.20 

- 0.50 

- 0.27 

- 0.22 

- 0.16 

- 1.00 

- 1.33 

- 0.80 

- 1.13 

- 1.61 

- 1.67 

- 1.42 

- 1.74 

- 0.17 

- 0.19 

- 0.08 

- 0.09 

- 0.39 

- 0.55 

- 0.39 

- 0.60 

- 0.65 

- 0.74 

- 0.76 

- 0.90 

- 0.23 

- 0.25 

- 0.10 

- 0.12 

- 0.50 

- 0.72 

- 0.51 

- 0,78 

- 0.85 

- 0.98 

- 0.98 

- 1 . 16 * 

- 0.15 

- 0.31 

- 0.08 

- 0.13 

1.46 

- 0.45 

- 0.41 

- 0.58 

0.00 

- 0.58 

- 0.67 

- 0 . 74 " 

- 0.21 

- 0.43 

- 0.11 

- 0.19 

1.46 

- 0.62 

- 0.57 

- 0.80 

0.00 

- 0.81 

- 0.92 

- 1.01 

- 0.10 

- 0.13 

- 0.06 

- 0.06 

- 0.21 

- 0.28 

- 0.22 

- 0.31 

- 0.32 

- 0.35 

- 0.37 

- 0.43 

- 0.12 

- 0.14 

- 0.06 

- 0.06 

- 0.22 

- 0.30 

- 0.22 

- 0.32 

- 0.34 

- 0.37 

- 0.37 

- 0.43 

- 0.19 

- 0.21 

- 0.10 

- 0.10 

- 0.37 

- 0.48 

- 0.37 

- 0.52 

- 0.56 

- 0.61 

- 0.62 

- 0.72 

- 0.27 

- 0.31 

- 0.14 

- 0.15 

- 0.53 

- 0.71 

- 0.54 

- 0.75 

- 0.82 

- 0.90 

- 0.89 

- 1.04 

- 0.39 

- 0.41 

- 0.19 

- 0.20 

- 0.75 

- 0.99 

- 0.74 

- 1.01 

- 1.14 

- 1.25 

- 1 . 2a 

- 1.40 

- 0.80 

- 0.91 

- 0.39 

- 0.42 

- 1.53 

- 2.02 

- 1.55 

- 2.12 

- 2.30 

- 2.53 

- 2.49 

- 2 . 88 * 

- 1.03 

- 1.15 

- 0.50 

- 0.56 

- 1.94 

- 2.57 

- 1.96 

- 2.68 

- 2.94 

- 3.23 

- 3.15 

-3.63 

- 1.14 

- 1.29 

- 0.57 

- 0.63 

- 2.18 

- 2.88 

- 2.20 

- 3.02 

- 3.30 

- 3.65 

- 3.56 

- 4.11 

- 0.10 

- 0.09 

- 0.04 

- 0.05 

- 0.21 

- 0.28 

- 0.20 

- 0.27 

- 0.31 

- 0.34 

- 0.32 

- 0.38 

- 0.10 

- 0.10 

- 0.05 

- 0.05 

- 0.21 

- 0.27 

- 0.19 

- 0.26 

- 0.31 

- 0.33 

- 0.32 

- 0.38 

- 0.28 

- 0.23 

- 0.13 

- 0.12 

- 0.54 

- 0.70 

- 0.50 

- 0.66 

- 0.82 

- 0.86 

- 0.82 

- 0.96 

- 0.43 

- 0.31 

- 0.20 

- 0.17 

- 0.81 

- 1.04 

- 0.71 

- 0.94 

- 1.21 

- 1.25 

- 1.15 

- 1.36 

- 0.70 

- 0.58 

- 0.32 

- 0.30 

- 1.31 

- 1.68 

- 1.20 

- 1.59 

- 1.95 

- 2.05 

- 1.91 

- 2.24 

- 0.96 

- 0.94 

- 0.46 

- 0.47 

- 1.81 

- 2.36 

- 1.73 

- 2.34 

- 2.70 

- 2.92 

- 2.78 

- 3.23 

- 1.15 

- 1.12 

- 0.55 

- 0.58 

- 2.17 

- 2.83 

- 2.07 

- 2.80 

- 3.25 

- 3.51 

-3.33 

-3.BT 

- 1.26 

- 1.26 

- 0.62 

- 0.64 

- 2.40 

- 3.13 

- 2.30 

- 3.13 

- 3.61 

- 3.91 

- 3.72 

- 4.33 

- 0.07 

- 0.05 

- 0.03 

- 0.02 

- 0.15 

- 0.20 

- 0.11 

- 0.18 

- 0.24 

- 0.25 

- 0.23 

- 0.28 

- 0.14 

- 0.09 

- 0.06 

- 0.05 

- 0.30 

- 0.41 

- 0.24 

- 0.36 

- 0.50 

- 0.51 

- 0.45 

- 0.55 

- 0.24 

- 0.15 

- 0.11 

- 0.08 

- 0.51 

- 0.70 

- 0.41 

- 0.61 

- 0.85 

- 0.88 

- 0.77 

- 0.96 

- 0.40 

- 0.22 

- 0.18 

- 0.13 

- 0.81 

- 1.09 

- 0.65 

- 0.92 

- 1.33 

- 1.37 

- 1.17 

- 1.43 

- 0.63 

- 0.39 

- 0.28 

- 0.21 

- 1.25 

- 1.67 

- 1.03 

- 1.44 

- 2.00 

- 2.08 

- 1.80 

— 2.13 

- 0.82 

- 0.62 

- 0.36 

- 0.32 

- 1.64 

- 2.20 

- 1.43 

- 2.00 

- 2.62 

- 2.75 

— 2 . 45 - 

- 2 . 94 ? 

- 0.98 

- 0.81 

- 0.44 

- 0.41 

- 1.99 

- 2.70 

- 1.76 

- 2.49 

- 3.22 

- 3.40 

- 3.06 

-3.ar 

- 1.07 

- 0.89 

- 0.48 

- 0.45 

- 2.18 

- 2.98 

- 1.93 

- 2.74 

- 3.56 

- 3.76 

-3.36- 

- 4.06 

- 0.02 

- 0.03 

- 0.01 

- 0.01 

- 0.06 

- 0.08 

- 0.06 

- 0.10 

- 0.09 

- 0.12 

- 0.13 

- 0.18 

- 0.03 

- 0.03 

- 0.01 

- 0.01 

- 0.07 

- 0.09 

- 0.07 

- 0.12 

- 0.11 

- 0.13 

- 0.14 

-o.ir 

- 0.06 

- 0.07 

- 0.03 

- 0.03 

- 0.14 

- 0.19 

- 0.14 

- 0.22 

- 0.23 

- 0.26 

- 0 . 28 * 

- 0.33 

- 0.14 

- 0.16 

- 0.07 

- 0.07 

- 0.32 

- 0.45 

- 0.32 

- 0.50 

- 0.53 

- 0.61 

— 0.63 

- 0.73 

- 0.21 

- 0.23 

- 0.09 

- 0.11 

- 0.47 

- 0.67 

- 0.48 

- 0.73 

- 0.78 

- 0.90 

- 0.91 

- 1.06 

- 0.29 

- 0.32 

- 0.13 

- 0.15 

- 0.63 

- 0.91 

- 0.65 

- 0.98 

- 1.07 

- 1.23 

- 1.23 

- 1.46 

- 0.76 

- 0.86 

- 0.34 

- 0.39 

- 1.56 

- 2.16 

- 1.58 

- 2.29 

- 2.50 

- 2.82 

- 2.77 

- 3.26 

- 0.85 

- 0.96 

- 0.39 

- 0.44 

- 1.75 

- 2.43 

- 1.77 

- 2.60 

- 2.83 

- 3.19 

- 3.16 

- 3.70 

- 0.03 

- 0.05 

- 0.02 

- 0.02 

1.46 

- 0.09 

- 0.08 

- 0.12 

0.00 

- 0.13 

- 0.13 

- 0.14 

- 0.03 

- 0.05 

- 0.02 

- 0.03 

1.45 

- 0.09 

- 0.08 

- 0.12 

0.00 

- 0.13 

- 0.13 

- 0.16 

- 0.05 

- 0.08 

- 0.03 

- 0.04 

1.46 

- 0.14 

- 0.12 

- 0.17 

0.00 

- 0.18 

- 0.19 

- 0 . 22 . 

- 0.12 

- 0.23 

- 0.07 

- 0.10 

1.46 

- 0.36 

- 0.32 

- 0.45 

0.01 

- 0.46 

- 0.53 

- 0.56 

- 0.19 

- 0.38 

- 0.10 

- 0.17 

1.46 

- 0.56 

- 0.51 

- 0.72 

0.00 

- 0.72 

- 0 . 82 , 

- 0 . 9 * 

- 0.31 

- 0.67 

- 0.16 

- 0.29 

1.46 

- 0.93 

- 0.86 

- 1.21 

0.00 

- 1.22 

- 1.37 

- 1.52 

- 0.72 

- 1.12 

- 0.36 

- 0.51 

1.46 

- 1.91 

- 1.62 

- 2.24 

0.00 

- 2.43 

- 2 . 56 - 

- 2.86 

V — 0.85 

- 1.27 

- 0.43 

- 0.59 

1.46 

- 2.21 

- 1 . B 7 

- 2.57 

0.00 

- 2.82 

- 2.94 

-3.31- 
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